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TIME-DEPENDENT  PROPAGATION  OF  HIGH-ENERGY 
LASER  BEAMS  THROUGH  THE  ATMOSPHERE:  II 

Abstract 

Various  factors  that  can  affect  coplanarity  should  benefit  multi- 

thermal  blooming  in  stagnation  zones  pulse  more  than  cw  beams.  The  methods 

are  examined,  including  stagnation- zone  of  treating  nonhorizontal  winds  hydro¬ 
motion,  longitudinal  air  motion  in  dynamically  for  cw  and  multipulse 

the  neighborhood  of  the  stagnation  steady-state  sources  are  discussed, 

zone,  and  the  effects  of  scenario  Pulse  Mself-b looming"  in  the  triangu- 

noncoplanarity .  Of  these  effects,  lar  pulse  approximation  is  discussed 

only  the  last  offers  any  reasonable  in  the  context  of  both  single  and 

hope  of  reducing  the  strong  thermal  multipulse  propagation.  It  is  shown 

blooming  that  normally  accompanies  that  self -blooming  and  multipulse 

stagnation  zones;  in  particular,  non-  blooming  cannot  be  treated  independently. 


1.  Introduction 


This  is  the  second  report  in  a 
series  dealing  with  the  general 
problem  of  time-dependent  thermal 
blooming  of  multipulse  and  cw  laser 
beams.'*'  Time  dependence  is  essential 
for  describing  the  propagation  of 
laser  beams  through  stagnation  zones, 
which  are  created  whenever  the  motion 
of  the  laser  platform  and  the  slewing 
of  the  laser  beam  combine  to  create 
a  null  effective  transverse  wind 
velocity  at  some  location  along  the 
propagation  path.  The  location  of 


vanishing  transverse  wind  we  shall 
call  the  stagnation  point,  and  the 
term  stagnation  zone  will  refer  to 
the  portion  of  the  propagation  path, 
extending  in  both  directions  from  the 
stagnation  point,  where  the  trans¬ 
verse  wind  has  not  yet  had  time  to 
blow  completely  across  the  beam.  The 
lack  of  wind  at  the  stagnation  point 
creates  a  steadily  decreasing  density 
and  a  thermal  lens  whose  strength 
grows  with  time.  This  report  will 
continue  the  study  of  stagnation 


zones  begun  in  Ref.  1,  and  discuss 
contributions  of  self-blooming  to 
multipulse  thermal  blooming  and  new 
models  that  have  been  added  to  the 
Four-D  code. 

Both  pivoted-absorption-cell 
2  3 

measurements 1 2 3 4  5  and  detailed  numer¬ 
ical  calculations  of  the  experimental 
1  4 

arrangements  9  give  evidence  that 

the  blooming  effects  of  stagnation 

zones  tend  to  saturate  with  time. 

Thus  the  beam  characteristics  seem 

to  approach  a  kind  of  quasi-steady 

state,  which  is  possibly  a  result  of 

the  steady  reduction  in  length  of  the 

3 

stagnation  zone  with  time.  Despite 
the  existence  of  these  quasi-steady 
states,  calculations  for  high-power 
beams  show  that  stagnation  zones  can 
lead  to  severe  beam  degradation. 

The  notion  of  a  stagnation  zone  re¬ 
quires  that  the  transverse  wind  veloc¬ 
ity  vanish  at  at  least  one  position 
along  the  propagation  path.  There 
are  always  present,  however,  a  number 
of  additional  effects  that  will  pre¬ 
vent  a  completely  stagnant  wind  con¬ 
dition  from  occurring  at  any 
position.  These  effects  are: 

1.  Natural  convection 

2.  Motion  of  the  stagnation  point 
with  time 

3.  Longitudinal  air  motion  at 
the  stagnation  point 

4.  Vertical  air  motion  due  to 

noncoplanar  scenario  geometry 


A  realistic  appraisal  of  the  influ¬ 
ence  of  stagnation  zones  on  beam 
propagation  requires  that  each  of 
these  effects  be  assessed  and  pos¬ 
sibly  incorporated  into  the  computa¬ 
tional  model. 

Natural  convection  flow  at  the 
stagnation  point  should  be  negligible 
for  practical  beam  sizes.  A  3.8-ym 
wavelength  and  a  474-kW  beam  power, 
for  example,  give  a  natural  con¬ 
vection  velocity  of  the  order  of 
10  cm/s.^  For  this  flow  velocity 
and  a  beam  radius  of  10  cm  at  the 
stagnation  point,  approximately  2  s 
is  required  for  the  beam  to  approach 
a  steady-state  density  distribution. 
This  time  is  excessive  for  preventing 
or  reducing  stagnation-zone  blooming 
effects,  which  may  develop  in  times 
ranging  from  1  ms  to  0.1  s.  At  a 
10.6-ym  wavelength  the  laser  heating 
rate  of  the  atmosphere  is  somewhat 
greater  for  the  same  intensity,  but 
the  natural  convection  velocity 
scales  only  with  the  cube  root  of  the 
absorbed  power,  so  flow  velocities 
would  not  be  significantly  above 
those  for  the  3.8-ym  case.  There¬ 
fore,  we  shall  not  consider  the 
effects  of  natural  convection 
further. 

Under  most  conditions  the  stag¬ 
nation  point  is  not  stationary  but 
moves  in  the  same  general  direction 
as  the  target  with  a  velocity  that 
is  little  different  from  the 


» 
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target's.  The  parcel  of  air  that 
sees  a  null  wind  speed  changes  with 
time  and  thus  does  not  heat  up  in 
the  manner  of  a  stationary  parcel. 

The  influence  of  this  stagnation- 
point  motion  on  beam  propagation  was 
found  to  be  minimal  for  a  cw  wave¬ 
form  example  treated  in  Ref.  1. 
Stagnation-point  motion  also  turns 
out  to  be  unimportant  for  a  multi¬ 
pulse  scenario  examined  in  this 
report.  The  conclusion  is  that 
stagnation-point  motion  is  unlikely 
to  have  much,  if  any,  effect  in 
alleviating  stagnation-zone  blooming, 
since,  despite  the  motion,  a  sub¬ 
stantial  propagation  path  exists  over 
which  wind  velocities  are  negligible. 

The  existence  of  a  null  transverse 
wind-velocity  component  at  a  stag¬ 
nation  point  in  no  way  guarantees  a 
vanishing  magnitude  of  the  wind  vec¬ 
tor,  because  a  nonvanishing  longi¬ 
tudinal  component  almost  always 
exists  there.  Any  air  parcel  found 
within  the  beam  at  the  stagnation 
point  will,  as  a  result,  exit  from 
the  beam  in  a  finite  length  of  time. 
Indeed,  in  coplanar  geometries  all 
wind-flow  trajectories  should  cross 
the  beam  in  two  locations:  one  for 
values  of  z  (longitudinal  position) 

below  the  stagnation  point  z  ,  and 

s 

the  other  for  values  above  z  .  The 

s 

wind  flow  may  be  in  either  the 
positive-  or  negative-3  direction. 

In  the  neighborhood  of  the  stagnation 


point,  the  wind-flow  trajectories 
will  enter  one  side  of  the  beam, 
reverse  direction  with  the  beam,  and 
exit  on  the  same  side.^  The  resi¬ 
dence  time  in  the  beam  for  fluid 
parcels  passing  through  the  beam 
center  at  the  stagnation  point  will, 
of  course,  depend  on  scenario  param¬ 
eters,  but  for  some  typical  beam 
sizes  and  scenarios  this  time  can  be 
of  the  order  of  0.5  to  1  s.  Longi¬ 
tudinal  flow  should  thus  be  as 
effective  as  natural  convection  in 
controlling  density  changes  at  the 
stagnation  point.  One  consequence 
of  these  re-entrant  wind-flow 
trajectories  is  that  air  densities 

for  z  values  greater  than  z  could 

s 

be  influenced  hydrodynamically  by 

densities  for  z  values  less  than  z  ; 

s 

but,  for  any  foreseeable  practical 
scenarios,  the  re-entrant  times  — 
except  perhaps  in  the  immediate 
neighborhood  of  the  stagnation 
point  —  would  be  considerably  longer 
than  times  of  interest.  Consequen¬ 
tly,  hydrodynamic  coupling  between 

points  above  and  below  z  can  be 

s 

safely  neglected  in  cases  of  prac¬ 
tical  interest. 

The  existence  of  a  position  where 
the  transverse  wind  velocity  vanishes 
presupposes  the  extremely  improbable 
coplanarity  of  the  laser  beam  and 
the  trajectories  of  its  platform  and 
the  target  —  a  situation  that  is 
clearly  a  limiting  case  of  real-world 
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scenarios,  which  are  invariably  non- 
coplanar.  In  the  more  general'  case  of 
noncoplanar  geometry,  only  the  wind 
component  along  a  certain  transverse 
axis  can  be  expected  to  vanish.  The 
wind  vector  in  the  transverse  plane 
will  rotate  and  attain  its  minimum 
magnitude  at  the  stagnation  point. 
Since  this  minimum  magnitude  can 
never  vanish,  except  in  a  space  of 
measure  zero,  a  steady  state  can 
always  be  defined  for  the  governing 
hydrodynamic  equations.  The  signif¬ 
icance  of  this  is  that,  in  systems 
analysis,  steady-state  numbers  can 
always  be  assigned  to  stagnation-zone 
situations,  at  least  for  some  nominal 
degree  of  noncoplanarity,  and  these 
numbers  can  be  obtained  from  simple 
and  relatively  cheap  steady-state 
calculations.  Truly  coplanar 
stagnation-zone  situations,  in  con¬ 
trast,  require  time-dependent  cal¬ 
culations  that  are  expensive  and 
require  considerable  care  in 
execution. 

In  the  coplanar  scenario  described 
in  Ref.  1,  for  example,  if  the  laser 
is  given  an  elevation  of  10  m  above 
the  plane  containing  the  target  and 
the  laser  platform,  the  vertical 
component  of  wind  velocity  at  the 
stagnation  point  takes  on  the  value 
of  1  m/s.  This  is  sufficient  to 
establish  a  steady  state  in  a  time 
of  the  order  of  0.1  s,  which  is  short 
compared  to  times  of  interest.  The 


small  vertical  velocity  component  at 
the  stagnation  point  leads  to  sub¬ 
stantial  changes  in  isointensity  con¬ 
tours  in  the  focal  plane,  but  the 
average  intensity  is  remarkably  close 
to  the  quasi-steady-state  value 
obtained  in  a  time-dependent  calcu¬ 
lation  for  the  corresponding  coplanar 
case. 

Thus,  small  amounts  of  non¬ 
coplanarity  should  not  be  expected 
to  greatly  improve  cw  laser  perfor¬ 
mance  in  stagnation-zone  situations 
but  should  contribute  to  ease  in 
understanding  and  predicting  it. 

The  case  of  multipulse  beams  is 
another  matter.  As  pulse-repetition 
frequencies  are  lowered,  a  small 
vertical  wind  component  at  the  stag¬ 
nation  point  becomes  more  and  more 
effective  in  sweeping  out  the  air 
between  pulses.  The  benefits  of  non¬ 
coplanarity  in  stagnation-zone 
situations  should  thus  be  greater  for 
multipulse  beams  than  for  cw  beams. 

The  current  status  of  the  Four-D 
code  is  summarized  in  Table  1,  and 
recent  additions  to  the  code  are 
described  in  the  body  of  this  report. 
We  have  continued  to  adhere  to  the 
philosophy  that  the  best  way  to 
approach  all  laser-propagation  cal¬ 
culations  is  through  a  single, 
unified  computer  code  that  can  be 
applied  to  any  problem.  The  advantages 
of  this  are  threefold.  First,  it 
greatly  simplifies  bookkeeping  (or. 
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more  appropriately,  code-keeping), 
since  a  proliferation  of  limited 
special-purpose  codes  is  avoided. 
Second,  if  each  type  of  calculation 
is  made  a  subset  of  a  larger  calcu- 
lational  capability,  new  features 
added  to  the  code  —  such  as  data- 


processing  routines,  adaptive-lens 
transformations ,  scenario  features , 
etc.  —  are  available  to  all  types  of 
calculation  at  once.  Third,  real¬ 
istic  simulations  are  possible,  since 
a  wide  range  of  conditions  can  be 
incorporated  into  any  calculation. 


Table  1.  Basic  outline  of  current  Four-D  propagation  code. 


Variables 


Form  of  propagation  equation 


a?*!/*** 

where  x3  y  are  transverse  coordinates  and 
z  is  axial  displacement. 

Scalar  wave  equation  in  parabolic  approxi¬ 
mation 


nk H  =  k2(-n2 "  d#  • 


Method  of  solving  propagation 
equation 


Hydrodynamics  for  steady-state 
cw  problems 


Transonic  slewing 


Symmetrized  split  operator,  finite  Fourier 
series,  fast  Fourier  transform  (FFT) 
algorithm 


x  =  k2(n2  -  1)  . 


Uses  exact  solution  to  linear  hydrodynamic 
equations.  Fourier  method  for  M  <  1. 
Characteristic  method  for  M  >  1.  Solves 


Steady-state  calculation  valid  for  all  Mach 
numbers  except  M  =  1.  Code  can  be  used 
arbitrarily  close  to  M  =  1. 
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Table  1  (continued) . 


Treatment  of  stagnation-zone 
problems  for  cw  beams 


Nonsteady  treatment  of  multi¬ 
pulse  density  changes 


Method  of  calculating  density 
change  for  individual  pulse 
in  train 


Treatment  of  steady-state 
multipulse  blooming 


Time-dependent  isobai;ic  approximation* 
Transient  succession  of  steady-state  den¬ 
sity  changes;  i.e.,  solves 


+  V 


Changes  in  density  from  previous  pulses  in 
train  are  calculated  with  isobaric  approxi¬ 
mation  using 

9pi  .  9pi 
dt  +  VX  dx 

~  2~ 01  ^  (•£>£/)  -  > 

Gs  n 

where  t In(x>y)  is  nth  pulse  fluence.  Den¬ 
sity  changes  resulting  from  the  same  pulse 
are  calculated  using  acoustic  equations  and 
triangular  pulse  shape. 

Takes  two-dimensional  Fourier  transform  of 


where  I  is  Fourier  transform  of  intensity, 
and  tp  is  the  time  duration  of  each  pulse. 
Source  aperture  should  be  softened  when 
using  this  code  provision. 

Previous  pulses  in  train  are  assumed  to  be 
periodic  replications  of  current  pulse. 
Solves 

3pl  _  3p1 

dt  +  Vx  dx  +  Vy  dy 


Treatment  of  turbulence 


i 


6  (t  -  t  ) 
n 


Pulse  self-blooming  is  treated  as  in  the 
nonsteady-state  case . 


Uses  phase-screen  method  of  Bradley  and 
Brown  with  Von  Karman  spectrum.  Phase 
screen  determined  by 


Table  1  (continued) . 


Lens  transformation  and 
treatment  of  lens  optics 


Treatment  of  nondiffraction 
limited  beams 


Adaptive  lens  transformation 


T(x,y) 


exp 


exp 


X 


a(k  ,k  )  §1/2(k  ,k  ) 
x9  yJ  n  x  x9  yJ 


where  a  is  a  complex  random  variable  and  $n 
is  spectral  density  of  index  fluctuations. 


Compensates  for  a  portion  of  lens  phase 
front  with  cylindrical  Talanov  lens  trans¬ 
formation,  Uses  in  spherical  case 


where  Zf  is  focal  length  of  lens,  Zf  is 
focal  length  compensated  for  by  Talanov 
transformation,  and  z £  is  focal  length  of 
initial  phase  front. 


Spherical-aberration  phase  determined  by 

.  SA  2t tA  ,  2  ,  2.2 

=  — o  (x  +  y  )  , 


or  phase-screen  method  of  Hogge  et  at* 
Phase  determined  as  in  turbulence,  only 


$  = 
n 


a2l2 
_ o 

2tt 


exp 


> 


2 

where  l0  is  correlation  length  and  a  phase 
variance. 


Removes  phase 
2 

X  [ai(*i  -  <xi>)2  +  h(xi  ~  <ai>)] 

1 

through  lens  transformation  and  deflection 
of  beam.  Here  x±  ~  x9  X2  =  2/,  averages  are 
intensity  weighted,  and  are  calcu¬ 
lated  to  keep  the  intensity  centroid  at 
mesh  center,  and  intensity  weighted  r.m.s. 
values  of  x  and  y  are  constant  with  s. 
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Table  1  (continued) . 


Selection  of  3-step 


Scenario  capability 


Treatment  of  multiline 
effects 


Treatment  of  beam  jitter 


Code  output 


Numerical  capacity  when  used 
with  CDC  7600  and  restricted 
to  internal  memory  (large 
and  small  core) 

Problem  zoning  features 


Adaptive  3-step  selection  based  on  limiting 
gradients  in  nonlinear  contribution  to 
phase.  Constant  3-step  over  any  portion  of 
range  also  possible. 

General  noncoplanar  scenario  geometry  capa¬ 
bility  involving  moving  laser  platform, 
moving  target,  and  arbitrary  wind  direction. 
In  coplanar  case,  wind  can  be  function  of 
t  and  3. 

Calculates  average  absorption  coefficient 
based  on  assumption  of  identical  field  dis¬ 
tributions  for  all  lines 

£  aifi exp  (_<V° 


exp(-0U2) 

i 

where  is  fraction  of  energy  in  line  %  at 
3  =  0. 

Takes  convolution  of  intensity  in  target 
plane  with  Gaussian  distribution: 


-  x\y  -  y')  dx'  d y'  , 

where  cr^  is  variance  introduced  by  jitter. 

Isointensity,  isodensity,  isophase,  and 
spectrum  contours.  Intensity  averaged  over 
contours.  Plots  of  intensity,  phase,  den¬ 
sity,  spatial  spectrum  along  specific 
directions,  etc.,  at  specific  times. 

Plots  of  peak  intensity  and  average  inten¬ 
sity  vs  time. 

Spatial  mesh,  64  x  64,  35  sampling  times, 
no  restriction  on  number  of  axial  space 
increments . 


Number  of  space  increments  in  x  and  y 
directions  must  be  equal  and  expressible  as 
a  power  of  2. 
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2.  Treatment  of  Moving  Stagnation  Zones 
in  Coplanar  Scenarios 


% 


The  basic  coplanar  scenario 
geometry  is  depicted  in  Fig.  1.  It 
is  assumed  that  the  target  and  laser 
platform  will  collide  on  the  x  axis 
after  a  time  Tq  has  elapsed.  The 
point  of  impact  is  denoted  by  P  and 
the  position  of  the  laser  by  L.  The 
effective  transverse  wind  speed  V, 
is  given  by  the  expression 


Vt(Te,a)  =  -sin  0y[7y 

-I  ('f  -  V  v»  Vi 

+  Vw  sinOj,  -  9y),  (1) 

where  V ^  is  the  target  (receiver) 
velocity,  is  the  laser-platform 
(transmitter)  velocity,  V ^  is  the 
background  wind  velocity,  R  is  the 


Fig.  1.  Diagram  of  coplanar  scenario  model. 
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range , 


Vp  ~  VE  “S  SS  •  <2> 

and 

Vn  -  VR  sin  0*.  (3) 


From  simple  trigonometry, 

..  f  cos  0  +  (7_t 

1  a  K  T  c 


Qt  =  cot 


sin  9 


a 


D  sin  d> 
_ _s 

sin  Qt  * 


] 


(4a) 

(4b) 


where  0  is  constant. 
a 

The  transverse  wind  speed  V^_ 

depends  on  and  hence  on  time 

through  the  dependence  of  the 

scenario  parameters  0^  and  R  on  time 

in  Eqs.  (4a)  and  (4b).  The  Four-D 

code  is  programmed  to  calculate 

VA T  j2)  as  a  function  of  T  and 
VO  o 

hence  of  time.  The  hydrodynamic 
equations  are  solved  numerically  by 
assuming  that  V ^  is  stepwise  constant 
over  each  integration  time  interval 
At. 

The  location  of  a  stagnation  point 
is  determined  by  setting  the  right- 
hand  side  of  Eq.  (1)  equal  to  zero. 


Clearly,  that  point  moves  with  time, 
and,  as  a  result,  a  different  parcel 
of  air  undergoes  heating  under  con¬ 
ditions  of  zero  wind  velocity  at  each 
instant.  In  general  the  stagnation 
point  will  move  with  a  velocity 
comparable  to  that  of  the  target. 

The  determination  of  from 
Eqs.  (l)-(4)  for  use  in  the  hydro- 
dynamic  equations  permits  an  accurate 
determination  of  the  effect  of  motion 
on  the  thermal  lens  in  the  stagnation 
zone.  It  is  more  difficult  to  follow 
the  irradiance  on  a  moving  target, 
however,  since  all  equations  are 
solved  in  a  retarded  time  frame.  It 
would  be  necessary  to  store  the 
irradiance  history  for  the  values  of 
z  corresponding  to  the  target  motion. 
Since  the  relative  change  of  the 
range  R  in  a  time  of  interest  is 
small,  it  is  a  good  approximation  to 
assume  the  focal  distance  and  the 
range  R  at  which  the  laser  intensity 
is  monitored  to  be  constants  in  time, 
while  the  correct  variable  R  is  taken 
into  account  in  the  hydrodynamic 
portion  of  the  calculation  by  means 
of  Eqs.  (l)-(4)  . 


« 


3.  Propagation  of  Multipulse  Laser  Beams 
Through  Stagnation  Zones 


The  particular  scenario  chosen 
leads  to  the  effective  transverse 
wind  as  a  function  of  range  for 


T  =  0,  shown  in  Fig.  2,  where  x  is 
measured  from  the  time  the  pulse 
turns  on.  The  pertinent  physical 
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data  are  the  following: 


Power,  P  53  kW 

Range,  R  2.5  km 

Focal  length,  /  4.5  km 

Absorption  coefficient,  a  0.25  km 

Wavelength,  X  10.6  ]im 

Aperture  diameter,  2 a 
2 

(Gaussian  at  1/e  )  30  cm 

Slewing  rate,  Q  7.44 

mrad/s 


Pulse-repetition  rate,  V  10,  25,  50, 

and  100  s”l 


These  data  can  be  expressed  in 
terms  of  the  following  dimensionless 
numbers : 


Np  =  kaZ/f  =  2.96, 

NQ  =  2a/vQht  =  3.2  x  10"2v, 
Ns  =  ttf/VQ  =  3.6, 

NA  =  af  =  1.125, 


where  N  ,  NOJ  NAJ  and  Nn  repre- 
r  O  b  A  U 

sent  respectively  the  Fresnel,  over¬ 
lap,  slewing,  absorption,  and  dis¬ 
tortion  numbers.  Although  the 
chosen  power,  53  kW,  is  rather  low, 
the  value  of  the  distortion  number 
is  quite  high,  and  the  resulting 
thermal  blooming  is  about  the  maxi¬ 
mum  that  the  code  can  accommodate. 

In  any  case,  the  above  parameters 


Axial  distance  —  km 

Fig.  2,  Transverse  wind  velocity  as  a 
function  of  axial  distance. 


are  adequate  for  assessing  the 
sensitivity  of  typical  multipulse 
laser  performance  to  stagnation-zone 
motion. 

In  the  present  scenario,  the 
stagnation  zone  and  target  move  with 
a  speed  of  300  m/s.  For  a  multi¬ 
pulse  beam  with  V  =  100  s  \  the 
stagnation  zone  moves  3  m  between 
pulses,  whereas  for  V  =  10  s  \  the 
stagnation  zone  moves  30  m  between 
pulses.  It  would  be  hoped  that,  in 
the  case  of  the  lower  pulse- 
repetition  frequency,  the  greater 
movement  of  the  stagnation  zone 
would  lead  to  a  reduced  buildup  of 
stagnant-air  density  changes .  This 
effect  turns  out  to  be  minimal. 

The  time  dependence  of  the  average 
intensity  on  target  (averaged  over 
the  minimum  half-power  area)  is 
shown  for  the  case  of  no  stagnation- 
zone  motion  in  Fig.  3(a)  and  for  the 
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Q  Vacuum,  corrected  for  linear  absorption 
a  Absorbing  atmosphere 


Area-averaged  target  intensity  as  a  function  of  pulse  time:  s: 


pulses  at  V  =  10  s-1.  (a)  No  stagnation-zone  motion  included 


(b)  Stagnation-zone  motion  included. 


0.4s 


0.5  s 


0.6  s 


Isointensity  contours  as  a  function  of  pulse  time  for  V  =  10  s 


Table  2. 

Comparison  of  multipulse  intensities  with  and 
zone  motion. 

without  stagnation- 

2 

Time-averaged  intensity  (W/cm  ) 

Time 

No  motion 

Motion 

(s) 

Average 

Peak 

Average 

Peak 

0.1 

168.5 

234.0 

168.5 

234.0 

0.2 

78.9 

117.3 

78.7 

117.3 

0.3 

60.7 

100.7 

60.8 

103.3 

0.4 

53.7 

92.5 

55.3 

96.8 

0.5 

49.6 

88.7 

52.4 

92.7 

0.6 

47.3 

84.1 

50.2 

87.5 

Maximum 

increase  resulting  from 

st agnation- zone  motion: 

av,  6.1%;  peak. 

4.0%. 

case  with  stagnation-zone  motion  in 
Fig.  3(b).  The  calculation  is  car¬ 
ried  out  for  six  pulses.  The  isoin¬ 
tensity  contours  for  the  six  pulses 
are  shown  in  Fig.  4  for  the  moving 
stagnation  zone.  The  contours  in 
the  nonmoving  case  are  so  similar 
that  they  are  not  shown.  The  per¬ 
formance  in  the  two  cases  is  sum¬ 
marized  pulse  by  pulse  in  Table  2, 
where  the  intensity  values  have  been 
averaged  over  the  interpulse 
separation  time,  and  the  percent 
improvements  in  intensity  indicated 
are  for  the  last  pulse  in  the  train. 

Surprisingly,  the  improvements 
in  peak  and  average  fluence  go  in 
the  opposite  direction.  The  peak 
and  average  fluences  are  actually 
slightly  higher  in  the  no-motion 
case,  as  shown  in  Table  3. 

This  behavior  is  due  to  the  large 
contribution  that  the  first  pulse 
in  the  train  makes  to  the  total 


fluence.  The  isofluence  contours 
for  the  case  with  motion  are  dis¬ 
played  in  Fig.  5.  The  central 
fluence  peak  contains  the  maximum 
value  and  makes  the  largest  contri¬ 
bution  to  the  fluence  averaged  over 
the  minimum  half-power  area. 
Apparently  in  the  no-motion  case 
the  subsequent  pulses  in  the  train 
make  a  greater  contribution  in  the 
central  region  than  do  the  corre¬ 
sponding  pulses  in  the  case  with 
motion.  This  small  difference  in 
peak  and  average  fluences  is  of 

Table  3.  Comparison  of  fluences  with 
and  without  stagnation-zone 


motion. 

Fluence 

(J/cm2) 

No  motion 

Motion 

Average  Peak 

Average  Peak 

30.61  42.1 

30.1  40.7 

Decrease  resulting 
zone  motion:  av,  1. 

from  stagnation- 
7%;  peak,  3.3%. 

-13- 


Fig.  5.  Fluence  contours  for  case  of 
V  =  10  s"1,  motion  included. 

little  or  no  practical  importance, 

and  is  indicative  of  the  fact  that 

stagnation- zone  motion  plays  no 

vital  role  in  determining  thermal 

blooming  in  stagnation  zones. 

Figure  6  shows  the  dependence  of 

average  intensity  on  time  at  the 

target  range  for  different  values 

of  pulse-repetition  rate,  V.  Each 

curve  begins  with  the  time  of 

arrival  of  the  second  pulse.  (The 

first  pulse  would  create  a  time- 

2 

averaged  intensity  of  189  W/cm  .) 

It  is  clear  that  reducing  V 


diminishes  the  effect  of  the  stag¬ 
nation  zone.  The  reason  obviously 
is  that  for  smaller  values  of  V  the 
air  can  be  swept  out  by  wind  between 
pulses  over  a  greater  proportion  of 
the  propagation  path. 

Sample  pulse-isointensity  con¬ 
tours  for  V  =  100  s  ^  are  displayed 
in  Fig.  7;  these  should  be  compared 
with  those  for  V  =  10  s  ^  in  Fig.  4. 
At  the  lower  repetition  rate,  the 
beam  has  divided  into  two  distinct 
spots.  At  the  higher  rate,  lateral 
peaks  are  also  formed  but  they  are 
much  less  distinct.  The  lateral 
spreading  of  the  contours  as  a 
function  of  time  is  shown  in  Fig.  8. 
The  width  perpendicular  to  the  wind 


Time  —  s 


Fig.  6.  Space-averaged  intensity  as  a 
function  of  pulse  time  for 
various  pulse-repetition 
rates . 
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,  7.  Isointensity  contours  for  pulses  sampled  from  V  =  100  s  train. 


0 


0.4 


Time  —  s 


Fig.  8.  Width  of  beam  in  direction 
perpendicular  to  the  wind 
at  various  pulse-repetition 
frequencies . 
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0.2 


0.6 


is  determined  by  measuring  the 
maximum  distance  perpendicular  to 
the  wind  direction  between  30%  con¬ 
tours  . 

In  conclusion,  the  performance 
of  a  multipulse  laser  under 


stagnation-zone  conditions  can  be 
improved  by  lowering  the  pulse- 
repetition  frequency,  but,  with  or 
without  motion  of  the  stagnation 
point,  the  thermal  blooming  is 
likely  to  be  substantial. 


4.  Effect  of  Longitudinal  Air  Motion  on  Flow 
in  the  Neighborhood  of  a  Stagnation  Zone 
for  Coplanar  Scenarios 


We  wish  to  examine  the  air- flow 
trajectories  in  a  coordinate  system 
that  moves  with  the  laser  beam. 

Take  the  x  axis  along  the  direction 
of  motion  of  the  laser  platform  and 
the  y  axis  perpendicular  to  it  in 
the  scenario  plane.  The  unit  vector 
z(t)  is  directed  along  the  rotating 
laser  beam,  and  £T(t)  is  taken 
normal  to  z(t)  (see  Fig.  9).  At 
any  instant  of  time,  z(t)  and  x1 (i) 
can  be  expressed  in  the  rest  frame 
of  the  laser  platform  by  means  of 
the  relations 

z(t)  =  [cos  0y(t) ,sin  0y( t)]9  (5a) 

x' (t)  =  [sin  0y(t) ,-cos  0y(£) ] > (5b) 
where 

0^  =  cos™*1  (B  •  x),  (6) 

and  the  angle  0^  is  calculated  from 
the  scenario  Eq.  (4a)  by  substi¬ 
tuting  T  -  t  for  T  .  The  common 

G  G 

origin  for  the  rest  frames  of  the 


laser  platform  and  the  rotating 
laser  beam  will  be  taken  at  the  cen¬ 
ter  of  the  laser  aperture. 

The  effective  wind  vector  in  the 
moving  coordinate  system  of  the 
laser  can  be  expressed  as 


L 


Fig.  9.  Vector  diagram  in  scenario 
plane;  z(t)  indicates 
instantaneous  direction  of 
slewing  laser  beam. 


Here  Vn  is  the  target  velocity  rela- 
— n 

tive  to  the  earth's  surface,  and 


^rel  =  4,  "  ’ 


(8) 


where  is  the  velocity  of  the 
laser  platform  and  is  the  veloc¬ 
ity  of  the  wind.  An  air  test  par¬ 
ticle  will  move  in  the  rest  frame 
of  the  laser  platform  along  trajec¬ 
tories  described  by 


r  (t)  =  r(0)  + 


(9) 


These  trajectories  can  be  expressed 
in  the  rotating  frame  of  the  laser 
beam  by  means  of  the  following 
relations : 


x'  it)  =  r(t)  •  x'  (t)  ,  (10a) 

z(t)  =  r_(t)  •  z(t)  .  (10b) 

Figure  10  shows  sample  air- 
particle  trajectories  in  the  vicin¬ 
ity  of  a  stagnation  point  for  four 
different  scenarios  of  practical 
interest  (of  the  type  shown  in 


Fig.  1).  The  target  range  and  stag¬ 
nation  point  location  at  t  =  0  are 
indicated  in  Table  4.  At  time  t  =  0 
the  particles  are  assumed  to  be 
located  precisely  at  the  stagnation 
point.  The  origin  of  the  transverse 
coordinate  is  assumed  to  be  at  the 
center  of  the  beam.  The  ticks  on 
the  trajectories  indicate  points 
separated  by  0.5  s  in  time.  The 
arrows  indicate  the  direction  of  air 
flow  with  increasing  time.  Also 
shown  in  Table  4  are  the  times  T 
actually  spent  in  the  laser  beam  by 
a  particle  that  crosses  the  stag¬ 
nation  point  at  the  center  of  a 
10-cm-radius  beam. 

The  longitudinal  wind  speed  in 

the  neighborhood  of  the  stagnation 

point  is  roughly  equal  to  as 

can  be  seen  from  Eq.  (7).  For  the 

scenarios  described  in  Table  4, 

V  ,  is  of  the  order  of  10  m/s.  For 
— : rel 

these  scenarios  the  longitudinal 
wind  component  will  be  of  limited  value 
in  clearing  the  beam  in  the  vicinity  of 
the  stagnation  point. 


Table  4.  Residence  time  in  beam  for  air  particles  passing  through  stagnation 
point. 


Scenario 

Target  position 
at  t  -  0 
(km) 

Stagnation  point 
at  t  =  0 
(km) 

Residence  time, 
T 

Cs) 

A 

1.5 

0.844 

1.0 

B 

1.0 

0.379 

0.5 

C 

2.5 

2.33 

1.7 

D 

1.0 

0.295 

1.0 
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Distance 


Transverse  coordinate  (x)  —  m 

Fig,  10.  Air-flow  trajectories  in  rest  frame  of  slewing  laser  beam  in  the 
presence  of  a  stagnation  zone.  Ticks  on  the  curves  denote 
0.5-s  intervals.  (a)  Range  =  1.5  km,  stagnation  point  at 
z  -  0.844  km.  (b)  Range  =  1.0  km,  stagnation  point  at 
z  =  0.379  km.  (c)  Range  =  2.5  km,  stagnation  point  at 

z  =  2.33  km.  (d)  Range  =  1  km,  stagnation  point  at  z  =  0.295  km. 


5.  Calculation  of  Transverse  Wind  Velocities 
for  Noncoplanar  Scenarios 


We  shall  again  assume  the 
scenario  of  Fig.  1,  only  now  we 
shall  relax  the* assumption  that  the 
scenario  or  kinematic  plane  neces¬ 
sarily  coincides  with  the  earth’s, 


or  the  horizontal,  plane.  The  line 
PLP  and  the  wind  vector,  however, 
will  be  assumed  to  lie  in  the  earth’s 
plane  (see  Fig.  11).  Again,  the  x 
direction  will  be  along  the  direction 
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of  motion  of  the  laser  platform,  the 
y  direction  will  be  in  the  kinematic 
plane,  and  the  unit  vector  normal 
to  the  kinematic  plane  will  be 


called  The  laser  aperture  will 

be  situated  at  position  L' ,  which 
is  at  a  height  h  above  the  line  PLP . 
The  line  LL'  defines  the  vector  li  = 
h  =  hh9  which  is  normal  to  the 
horizontal  plane  and  makes  an 


angle  0  with  the  vector  The 

r  __ 

scenario  parameters  D >  R* 

and  are  now  defined  in  a  plane 
tilted  with  respect  to  the  horizontal 
plane,  but  they  are  related  exactly 
the  same  as  before.  The  distance  i?, 
nowever,  no  longer  has  the  signifi¬ 
cance  of  range.  The  calculation  of 


the  true  range  R '  is  described  below. 

In  order  to  follow  the  wind  in  a 
frame  of  reference  that  moves  with 
the  laser,  it  is  necessary  to  intro¬ 
duce  an  appropriate  orthogonal 
coordinate  system.  Clearly  this 
coordinate  system  will  not  be 
unique,  but  a  suitable  one  can  be 
defined  as  follows:  z  is  directed 
along  the  laser  beam, 


and 


y 


x  v 
-T 

X  IT\’ 


(11a) 


x'  =  x  g  .  (11b) 


It  is  most  convenient  to  express  all 
vectors  used  in  the  computation  in 
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9 


(12f) 


the  kinematic  coordinate  system. 


Hence 

we  have 

^  > 

1! 

(1,  0,  0)  , 

(12a) 

II 

< 

(cos  9^,  sin  0^,  0), 

(12b) 

t=Q  > 

II 

(cos  9y,  sin  9y,  0), 

(12c) 

✓s. 

h  = 

(0,  sin  0^,  cos  0^), 

(12d) 

ii 

< 

(cos  0^,  sin  9^  cos  0^, 

-  sin  0r.  sin  0  ) 

W  p 

(12e) 

R'  =  R  -  h 

/\  /\  /\ 

where  V  ^ W  are  un*t  vectors 
directed  along  Vy,  V^,  and  V^,  and 
the  vectors  R  =  RR  and  R'  are  directed 
along  lines  extending  from  L  and  L ' , 
respectively,  to  the  receiver 
(target) . 

The  effective  wind  seen  in  the 
frame  of  reference  moving  with  the 
laser  beam  is 


veff(s,t)  =  (v,,  -  V^) 


■^w  - T ' 


-f{' 


-  [v^  - 


v-4? 


>}■ 


(13) 


The  effective  horizontal  and  ver¬ 
tical  wind  components  in  Eqs.  (14) 
become  inputs  to  the  hydrodynamic 
calculation  which  is  described  in 
Sections  6  and  7. 

6.  Steady-State  Solutions  of  Hydrodynamic 
Equations  for  Arbitrary  Transverse  Wind  Velocities: 

cw  Steady  State 

Noncoplanar  scenarios  create 
effective  winds  whose  orientation 
in  the  transverse  plane  vary  with 
propagation  distance  a.  All 
symmetry  in  the  transverse  plane 
is  lost,  and  the  x  axis  can  no 
longer  serve  as  the  wind  axis. 

The  linearized  hydrodynamic 
equations  must  be  recast  and 
solved  for  a  wind  having  an 
arbitrary  direction. 


The  linearized  hydrodynamic 
equations  to  be  solved  are 

dPi 


d  t 


p0  -  2.1  =  °> 


(15a) 


p0  dt  -1  ?Pi  +  to 


f 3y 


n 


1  i 


to 


,  3ylt7  2  «  v 
+  ~  — —  —  o  .  .V  •  v, 


9a; . 


3  13—  -1/ 


(15b) 


d  ,  2  . 

dt  (pl  ~  %  pl}  =  (Y  ~  1)  al,  (15c) 


The  effective  wind  components  along 
x’  and  y '  are  then  obtained  from 


^effm'  -^eff 


V  «  ,  =  V  „ 
effy'  — eff 


r 


(14a) 

(14b) 
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where  p^,  v^,  and  p^  represent  the 
density,  velocity,  and  pressure  per¬ 
turbations  induced  by  laser  heating, 
rt  is  the  viscosity,  and  the  total 

derivative  77  is  defined  by 
at 


d  3  ,  3  ,  3  r\ 

-jT  =  ttt  +  V  7 - h  y  7 — .  (16) 

d t  3 t  x  dx  y  3z/ 


Elimination  of  p^  and  yields  the 


following  equation  for  : 


d  [  d^  ^  „  2^2  4  n  ry2 

d t  \  dt2  P1  "  s  V  P1  “  3  pQ  V  P1J 


=  (Y  -  1)  av  I  .  (17) 

We  are  interested  in  tHe  steady  state 
or  the  case  in  which  Eq.  (17)  becomes 


.  9  .  3 

Vx  dx  Vy  3  y 


■  2W2 
-  V 


1  pj  ^(V,h+”yJyi 


Pi  =  (Y  -  l)otV  I,  (18) 


where 


3Pn 

P1  Vx  3a; 


3pn 


+  v 


y 


(19) 


The  solution  for  p^  is  carried  out  in 

two  steps:  first  Eq.  (18)  is  solved 

with  p^  as  dependent  variable,  and 

then  Eq.  (19)  is  solved  for  p^. 

We  shall  restrict  our  attention  to 

2  2 

the  subsonic  case,  where  V  +  V 

2  *  y . 

<  o  .  In  that  case  Eq.  (18)  is  ellip- 
s 


tic  and  can  be  expressed  in  terms  of 
a  finite  Fourier  series  representation 


P-,  C x»y ) 


=  S  pi(vy  exp[^(y +  y)]- 

(20) 


k  ,k 
x  y 


The  coefficients  p(k^k  )  satisfy 


\{k  ,k  )  =  - 

lv  x’  y ' 


aT(kx,ky)(kl  4 


(k^  +  k b  1  +  4-  £  p9  -(p  k  +  v  k  )  - 

«  z/'L  3  p0e2  a;  a;  y  y' _ 


(v  k  +  v  k  V 

1  L-ud. 

2 


(21) 


'V 

where  _T(fc  tk  )  is  the  Fourier  trans- 
x  y 

form  of  I(x,y).  The  inverse  Fourier 
transform  of  Eq.  (21)  is  evaluated  by 


means  of  the  fast  Fourier  transform 
(FFT)  algorithm.  The  function 
p ^(x,y)  thus  obtained  then  becomes 
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If  we  define 


the  source  term  for  Eq.  (19),  which 
can  be  solved  using  the  Carlson  method"*- 
for  integration  along  characteristics. 
The  solution  of  Eq.  (19)  is  obtained 
by  a  difference  method  in  configura¬ 
tion  space  in  preference  to  a  Fourier 
transform  method  because  the  transform 
of  p ^(x,y)  will  have  poles  whenever 
Vx^x  +  Vy^-y  =  0*  The  evaluation  of 
the  inverse  transform  by  the  FFT 


V 

JL 

A 

X 

V 

X 

A 

y 

V 

X 

9 

V 

_JL 

5 

'V 

(22a) 


(22b) 


(22c) 


algorithm  will  be  troublesome,  since  the  difference  equation  satisfied 


these  poles  must  be  avoided. 


by  -  p^(iAar, jAz/)  can  be  written 


J  •  •  —  -jr  w  • 

id  3  t 


L~i '  >  d~d  '  +  (X  3)  PA,J-J  ' 

uT\  +  i  +  (x "  i)  pa,j-j'}  for 


3  >  1; 


(23a) 


PAj  B  Pi-i\d~d'  +  (1  Pi-V  J 


+  2\vx\  yid  +  B  Pi-i',d-d'  +  (1  ‘  B)  Pi-i ' ,  j  j  for  6  -  1* 


(23b) 


7.  Steady-State  Solutions  of  Hydrodynamic  Equations  for 
Arbitrary  Transverse  Wind  Velocities: 
Multipulse  Steady  State 


Isobaric  density  changes  induced  by  multipulse  heating  are  governed  by 
the  equation 


3e”P  8p"p  3p”I 

+  VX  dx  Vl j  ty 


“  01  S  TpJn(x'y)  6(t  -  V  5 


where  represents  the  fluence  of  the  nth  pulse,  and  x  represents 
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the  pulse  width.  If  "steady-state" 
conditions  prevail,  it  can  be  assumed 
that  does  not  vary  from  pulse  to 


pulse.  Hence  I  (x9y)  =  I(x3y) . 

Taking  the  Fourier  transform  of  Eq.  (24) 
with  respect  to  x  and  y  yields 


3p 


mp 


'w  +  i!-vxkx  +  = "  3L=r’  “V(VV  (25) 

Q 


n 


Solving  Eq.  (25)  for  at  a  time  t  =  mht,  where  m  is  any  integer  and  At 
is  the  time  interval  between  successive  pulses,  gives 


N 


Y  -  1 
2 


L 

a  T(k  ,k  )  Y  exv[-inAt(k  v  +  k  v  )] .  (26) 

x  y  y  y 


n- 1 


The  exponentials  in  Eq.  (26)  corres¬ 
pond  to  translations  of  the  individual- 
pulse  fluence  distributions  in  con¬ 
figuration  space  by  wind  motion. 

The  summation  begins  with  n  =  1 
because  the  isobaric  density  changes 
created  by  a  given  pulse  do  not  have 
time  to  develop  during  the  pulse 
width,  T^.  The  upper  limit  N  is 
based  on  numerical  considerations  and 
is  determined  by 


[-#d  -  ™ 

where  N  is  the  numerical  length  of  the 
mesh  used  for  solving  the  wave 
equation. 


In  Eq.  (27)  MIN  signifies  the  min¬ 
imum  of  the  arguments,  and  the  square 
brackets  represent  the  integer  part 
of  the  arguments  inside  them.  An 


stagnation  point  is  encountered  along 
the  propagation  path.  In  such  cases 
the  total  density  change  at  the  stag¬ 
nation  point  can  be  kept  bounded. 

For  example,  ^npUt  might  be  set  equal 
to  the  actual  number  of  pulses  in  a 
given  train,  in  which  case  a  true 
lower  bound  could  be  assigned  to  the 
intensities  at  the  target. 

The  remaining  arguments  in  Eq.  (27) 
prevent  any  pulse  fluence  distribution 
from  affecting  the  density  calculation 
if  it  has  been  translated  by  more 
than  the  minimum  (physical)  dimension 
of  the  computational  mesh  for  the 
wave  equation,  i.e.  MIN (Mar,  My). 


-23- 


The  density  calculation  itself  is  of  the  density  contributions  by  past 

carried  out  on  a  2 N  x  2 N  mesh-,  which  pulses  in  the  train  is  avoided, 

has  a  buffer  of  length  N  in  both  the 

x  and  y  directions.  Thus  if  Np  The  summation  in  Eq.  (26)  may  be 

satisfies  condition  (27),  periodic  evaluated  directly,  and  p"  can  be 

"wrap-around"  or  positional  aliasing  expressed  in  the  form 


The  density  p^(x,y)  is  then  obtainable  may  not  correspond  to  lattice 

from  Eq.  (28)  by  an  inverse  transform  translation  operators  on  the  compu- 

operation  using  the  FFT  algorithm.  tational  mesh.  In  such  a  case 

Equation  (28)  has  been  used  in  a  num-  ringing  can  be  suppressed  by  express- 

ber  of  test  examples  with  satisfactory  ing  the  solution  of  Eq.  (26)  in 

results.  If  the  spectrum  I(k.k)  is  terms  of  the  interpolations  of  lattice 

jj  y 

particularly  rich  in  high  spatial  fre-  shift  operations. 


quencies,  Eq.  (28)  may  give  rise  to  a 
ringing  behavior  in  configuration  space 
due  to  the  fact  that  the  shift  operators 

exp  (-inktkv),  exp  {-inhtk  v  )  (29) 

JC  X  y  y 


By  means  of  bilinear  interpolation, 
one  can  express  any  function  T(x^y) 
at  positions  intermediate  to  the  lat¬ 
tice  by  means  of 


rj+fx>k+fy  <i  4*  4  Tc+i,k +  a  yj  fy  Tj,k+ 1 


+  44  4 


y  Tc*iMi  +  a  '  4K1 ' 4>  Tj,k  ■  <30> 


0  <  4  1 1  .  0  <  fy  <  1  , 
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where  f  and  f  represent  fractional  may  thus  be  represented  in  the  fol- 
x  y 

distances  between  lattice  coordinates,  lowing  alternative  form,  which  avoids 
and  where  the  numbers  T .  ,  represent  the  use  of  nonlattice  shift  operators 


values  of  T(x^y)  sampled  at  lattice 
points.  The  summation  in  Eq.  (26) 


(the  notation  [  ]  signifies  the 
integer  part  of  the  argument) : 


N 


2  exp  [~inht{kxvx  +  kyVy)] 


n- 1 


N  r 
P  r 


-  1 


n=l 


4<X 


fy )  exp(i  ^ 


([n  n]  +  1)  +  n.  [Tin] 


X 


’a?""  '  '  k  -  Y 

2/  * 


-» 


+  V1  _  v  exp  v  £  {\[v]  +  %([y]  + 1} 


■}) 


+  Vy  exp 


f  *  {" 


k  ([ry*]  +  1)  +  ( [n  n]  +  1) 

a?  ^ 


}) 


+  (1  -  (1  ~  jp  exp  ( y 


£  I  {"* 1 V1  +  nk  tv1}) 


r  -  ?<w  <3« 


where 


y  At 
ar 

'a;  Aa: 


- 


y  At 

ny  =  * 


4(n)  =  v  -  ty*1* 


4(n)  =  y "  [y]  • 


The  summation  in  Eq.  (26)  over  n  can 
be  evaluated  as  a  2/17  x  2/17  DFT  with 
the  aid  of  the  FFT  algorithm.  For  a 
given  value  of  w,  the  numbers  [rj  n] , 

X 

[t\  n]  +  1  can  each  be  identified  as 

SC 

^-coordinates  n  and  the  numbers 
x 

(32)  y >  fn^w]  +  1  as  ^/-coordinates  n 
in  the  lattice  space.  Thus  each 
exponential  in  the  summation  in 
Eq.  (31)  can  be  identified  with  a  par¬ 
ticular  lattice  point  n  ,  n  .  As  the 

x  y 

index  n  is  incremented,  the  appro¬ 
priate  bilinear  function  of  f  (n)  and 

x 
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f  (w)  is  added  to  the  contents  of  a 

y 

storage  register  corresponding  to 

coordinates  n  ,  n  .  On  completion  of 

cc  z/ 

this  operation,  a  two-dimensional  DFT 
of  the  resultant  array  will  yield  the 
desired  sum  (31). 

The  Fourier  transform  of  the 
density  can  then  be  expressed  as 

\{kx,kyimht)  = 

a 

x  «?<yy  ?<vy  •  <33) 


Both  the  options  (28)  and  (29)  are 
currently  available  in  the  Four-D 
code,  and  the  cases  run  have  produced 
results  that  are  almost  indistinguish¬ 
able. 

The  shifts  and  interpolations 
implied  in  Eq.  (31)  may,  of  course, 
be  carried  out  strictly  in  configura¬ 
tion  space.  If  N  is  small,  this 

tr 

procedure  may  be  more  economical. 

As  N  becomes  large,  the  Fourier 
transform  method  becomes  more 
economical. 


8.  Effect  of  Noncoplanarity  on  Propagation  of 
cw  Laser  Beams  Through  Stagnation  Zones 


We  shall  focus  attention  on  the 
scenario  discussed  in  Ref.  1,  in 
which  the  total  propagation  dis¬ 
tance  is  1.5  km  and  the  stagnation 
point  occurs  at  z  =  0.8439  km.  The 

initial  diffraction-limited  beam  is 
2 

Gaussian,  with  1/ e  -intensity 
diameter  of  70  cm,  and  is  assumed 
to  be  focused  at  the  1.5-km  range. 

The  wavelength  and  absorption  coef¬ 
ficient  are  assumed  to  be  3.8  ym 
and  0.07  km  ^ ,  respectively.  For 
reference  the  results  of  the  time 
dependent  calculations  at  t  =  60  ms 
are  given  in  Table  5.  For  this 
value  of  t9  the  beam  properties  are 
changing  very  slowly,  and  the  assump¬ 


tion  of  a  "quasi"  steady  state  is  a 
reasonable  one. 

In  the  noncoplanar  scenario,  on 
the  other  hand,  a  true  steady  state 
is  known  to  exist,  and  a  time  to 
establish  this  steady  state  can  be 
estimated  by  dividing  the  beam 
diameter  by  the  magnitude  of  the 
vertical  wind  component  at  the 
stagnation  point.  The  noncoplanar 
results  are  naturally  much  cheaper 
to  obtain  than  the  corresponding 
coplanar  results. 

In  Table  6,  steady-state  results 
are  given  for  the  scenario  corres¬ 
ponding  to  Table  5  for  a  variety  of 
elevations  of  the  laser  aperture 
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Table  5.  Beam  properties  on  target  at  t  =  60  ms. 


Laser  power 
(kW) 

Peak  intensity 
at  target 
(kW/cm^) 

Minimum  half¬ 
power  area 

(cm^) 

Intensity  averaged 
over  minimum 
half-power  area 
(kW/cm^) 

500 

10.8 

33.4 

6.72 

500a 

9.8 

53.5 

4.19 

500b 

11.0 

33.2 

6.76 

250 

12.4 

13.5 

8.34 

125 

17.7 

4.42 

12.7 

62.5 

22.7 

1.65 

17.0 

aFocus  100 
^Motion  of 

m  beyond  range. 

stagnation  zone  taken 

into  account. 

Table 

6.  Steady-state  cw  beam  properties  as  a  function  of 
scenario  plane. 

laser  height  above 

Laser  Laser 

power  elevation 
(kW)  (m) 

Vertical  wind 
speed  at 

stagnation  point 
(m/s) 

Minimum  half- 
power  area 
(stagnation 
point) 

(cm^) 

Minimum  half¬ 
power  area 
(target) 

(cm^) 

Time  to 
steady 
state 
(s) 

Peak 

intensity 
at  target 
(kW/cm^) 

Intensity 
averaged  over 
minimum  half¬ 
power  area  at 
target 
(kW/cm^) 

500 

5 

0.55 

293 

37.7 

0.312 

11.0 

5.97 

10 

1.1 

291 

33.6 

.155 

12.0 

6.69 

20 

2.2 

290 

29.1 

.077 

14.0 

7.72 

30 

3.3 

289 

26.4 

.052 

15.9 

8.52 

40 

4.4 

287 

23.6 

.039 

16.1 

9.53 

250 

5 

0.55 

279 

13.7 

.303 

16.6 

8.2 

10 

1.1 

278 

12.0 

.152 

17.8 

9.36 

20 

2.2 

277 

10.1 

.076 

19.5 

11.0 

30 

3.3 

278 

8.9 

.051 

21.6 

12.6 

40 

4.4 

276 

7.88 

.038 

24.1 

14.3 

125 

5 

0.55 

272 

4.66 

.300 

22.9 

12.0 

10 

1.1 

271 

4.04 

.150 

25.3 

13.9 

20 

2.2 

271 

3.33 

.075 

28.6 

19.5 

40 

4.4 

270 

2.31 

.037 

39.7 

24.3 

62.5 

5 

0.55 

268 

1.58 

.29 

30.6 

17.7 

10 

1.1 

268 

1.37 

0.14 

34.3 

20.5 
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Magnitude  m/s  y  component  m/s  x  component 


.  C  _  Stagnation 
point 


0  0.5 


0  0.5  1.0  1.5 


0  0.5  1.0  1.5 

Axial  distance  —  km 

Fig.  12.  Transverse  wind  velocity  as 
a  function  of  axial  distance 
for  cw  beam.  (a)  x  com¬ 
ponent.  (b)  y  component. 

(c)  Magnitude. 


above  the  scenario  plane.  Figure  12 
shows  the  variations  with  z  of  the 
horizontal  and  vertical  components 
and  the  magnitude  of  the  transverse 
wind. 

From  Tables  5  and  6,  it  is  evi¬ 
dent  that  the  space-averaged  inten¬ 
sities  in  the  focal  plane  for  the 
noncoplanar  scenario  at  5-m  eleva¬ 
tion  agree  with  the  corresponding 
average  steady-state  intensities  for 


U>0 


M 

\ijf\f- l  \  l 

mm 


W  W  ^  i  M\ 

II  •Iff 

n# 


Time -depen dent 
coplanar,  quasi¬ 
steady  state 


Steady  state, 
noncoplanar, 
h  =  10  m 


Fig,  13,  Comparison  of  isointensity 
contours  for  stagnation- 
zone  situations  in  coplanar 
and  noncoplanar  cases. 


the  coplanar  scenarios  to  within 
less  than  10%.  The  peak  intensities 
for  the  noncoplanar  scenario  at  5  m, 
on  the  other  hand,  are  somewhat 
higher  than  the  corresponding  values 
for  the  coplanar  case.  There  is 
also  a  substantial  difference  in  the 
appearance  of  the  isointensity  con¬ 
tours  in  the  focal  plane  (Fig,  13). 
As  would  be  expected,  performance 
improves  with  height,  although  the 
improvement  is  marginal  for  the 
elevations  considered.  In  all  cases 


a  steady-state  condition  can  be 
reached  in  a  time  small  compared  with 
times  of  interest. 

In  conclusion,  average  intensi¬ 
ties  for  coplanar  stagnation-zone 
scenarios  can  be  calculated  by 
adding  nominal  noncoplanar  features 
to  the  scenario  and  performing  a 
steady-state  calculation.  For 
cw  beams,  however,  rather 
substantial  laser  elevations 
must  be  provided  to  alleviate 
stagnation-zone  effects. 


9.  Effect  of  Noncoplanarity  of  Propagation  of 
Multipulse  Beams  Through  Stagnation  Zones 


We  turn  our  attention  again  to 
the  scenario  of  Section  3.  All 
problem  parameters  are  the  same, 
except  that  the  laser  is  now  assumed 
to  be  elevated  10  m  above  the 
scenario  plane.  Figure  14  shows  the 
vertical  and  horizontal  components 
of  transverse  wind  velocity  as 
functions  of  propagation  distance. 
Figure  15  shows  the  isointensity 
contours  in  the  target  plane  for  the 
various  repetition  rates. 

Table  7  compares  laser  perfor¬ 
mance  on  target  as  a  function  of 
pulse-repetition  frequency  for  the 
coplanar  scenario  and  the  noncoplanar 
scenario  with  a  laser  elevation  of 
10  m.  In  the  absence  of  complete 


steady-state  data  for  the  coplanar 
case,  we  have  used  in  Table  7  inten¬ 
sity  values  corresponding  to  the 
final  times  exhibited  in  Fig,  6  for 
a  given  value  of  V.  Thus  the 
improvements  due  to  noncoplanarity 
shown  in  Table  7  are  conservative 
estimates . 

It  is  seen  from  Table  7  that 
improvements  of  at  least  a  factor 
of  2,  conservatively  estimated,  are 
possible  for  all  values  of  V.  In 
the  case  of  V  =  10  s~^  the  laser 
performance  is  even  better  than  it 

would  be  in  a  vacuum.  The  reason 
is  that  for  this  pulse-repetition 
frequency  the  overlap  number  at  the 
stagnation  point  is  only  2,  and  for 
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Axial  distance  —  km 

Fig.  14.  Transverse  wind  velocity  as 
a  function  of  axial  distance 
for  multipulse  beam.  (a)  x 
component.  (b)  y  component. 


overlap  numbers  in  the  range  1-2 
such  enhancement  effects  for  multi- 

g 

pulse  beams  are  well  known. 


Fig.  15.  Changing  shapes  of  isoin¬ 
tensity  contours  as  a 
function  of  pulse-repetition 
rate  for  noncoplanar  scenario 
laser  at  10-m  elevation. 


To  summarize:  there  is  clearly 
some  hope  of  minimizing  stagnation- 
zone  blooming  for  multipulse  beams 
by  a  combination  of  elevating  the 
laser  aperture  above  the  scenario 
plane  and  lowering  the  pulse- 
repetition  frequency. 


[i] 


Table  7.  Comparison  of  multipluse  beam  properties  for  coplanar  and  non- 

coplanar  scenarios.  Power  =  53  kW,  range  =  2.5  km,  A  =  10.6  pm, 
elevation  h  -  10  m,  and  vertical  wind  speed  at  stagnation  point 
=  0.61  m/s. 


Pulse 

repetition 

frequency, 

V 

(s-1) 

Minimum 

half-power 

area 

(stagnation 
point , 
noncoplanar 
scenario) 
(cm”2) 

Time  to 
steady 
state 
(non¬ 
coplanar 
scenario) 
(s) 

Overlap 
number  at 
stagnation 
point  (non¬ 
coplanar 
scenario) 

Peak 

intensity 
at  target 
(coplanar 
scenario) 
(W/cm2) 

Peak 

intensity 
at  target 
(non¬ 
coplanar 
scenario) 
(W/cm2) 

Intensity 

averaged 

over 

minimum 
ha If -power 
area 

(coplanar 

scenario) 

(W/cm2) 

Intensity 

averaged 

over 

minimum 
half -power 
area  (non¬ 
coplanar 
scenario) 
(W/cm2) 

10 

131 

0.19 

1.9 

85. 5C 

287a 

52.0 C 

181b 

25 

116 

0.18 

4.49 

59. 5C 

116 

32. 5C 

65.6 

50 

104 

0.17 

8.49 

28. 7d 

70.9 

17. 8d 

42.3 

100 

140 

0.19 

19.0 

30. 4d 

49.0 

13. 2e 

30.3 

Vacuum  beam  has  value  238* 


^Vacuum  beam  has  value  170. 

C£  -  0.6  s,  steady  state  has  not  been  reached. 

=  0.32  s,  steady  state  has  not  been  reached. 
=  0.2  s,  steady  state  has  not  been  reached. 


10.  Single-Pulse  Thermal  Blooming  in  the 
Triangular  Pulse  Approximation 


The  isobaric  approximation  for 
changes  in  air  density  is  invalid 
for  a  single  laser  pulse  whose  dura¬ 
tion  is  comparable  to  or  less  than 
the  transit  time  of  sound  across  the 

beam.  In  this  time  regime  —  referred 
3 

to  as  the  t  -regime  because  of  the 
time  dependence  of  density  changes 
arising  from  an  applied  constant 
laser-energy  absorption  rate  —  the 
air-density  changes  must  be  deter¬ 
mined  from  the  complete  set  of  time- 
dependent  hydrodynamic  equations, 

Eqs .  (15  ^ 


At  late  times  in  the  pulse,  t 
thermal  blooming  tends  to  reduce  the 
on-axis  intensity  relative  to  what 
it  would  be  if  the  beam  were  propa¬ 
gating  in  vacuum.  This  reduction 
increases  with  time,  and  for  suf¬ 
ficiently  late  times  a  depression 
appears  in  the  center  of  the  beam. 
Energy  added  to  the  pulse  at  later 
times  will  contribute  only  margin¬ 
ally  to  the  on- axis  fluence.  Thus, 
for  a  specific  peak  pulse 
intensity,  the  on- axis  fluence 
appears  to  saturate  as  the  pulse 
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duration  is  stretched  out  more 
and  more. 

These  properties  are  best  illus¬ 
trated  by  a  numerical  example.  Let 

us  consider  a  beam  that  is  Gaussian 
2 

at  z  =  0  with  1/e  -intensity  radius 
25  cm.  The  beam,  which  is  focused 
at  2.5  km,  is  assumed  to  be  2x  dif¬ 
fraction  limited  (A-scaled)  with 
X  =  10.59  ym  and  a  =  0.3  x  10  ^  cm  \ 
The  pulse  is  square-shaped  in  time 
and  lasts  100  ys.  The  choice  of  a 
square-shaped  pulse  is  convenient 
because  a  single  calculation  con¬ 
tains  the  complete  information  for 
all  square  pulses  of  duration 
shorter  than  the  one  chosen. 

Figure  16  shows  the  on-axis 
intensity  at  2  =  2.0  km,  obtained 


Fig.  16.  On-axis  intensity  as  a 

function  of  time.  The  pulse 
is  taken  to  be  square-shaped 
in  time.  Thermal  blooming 
reduces  on-axis  intensity 
to  a  negligible  value  after 
a  sufficiently  long  time. 


by  detailed  numerical  solution^  of 
Eqs.  (15).  The  on-axis  intensity 
clearly  drops  to  a  negligible  value 
before  the  end  of  the  pulse,  and, 
as  a  consequence,  the  on- axis 
fluence  saturates  as  the  pulse  width 
increases,  as  can  be  seen  in  Fig.  17. 
The  detailed  temporal  evolution  of 
the  spatial  shape  of  the  beam  is 
shown  in  Figs.  18  and  19.  Figure  18 
is  a  three-dimensional  plot  of  the 
laser  intensity  as  a  function  of 
time  and  radius.  Figure  19  shows 
the  radial  intensity  profiles  for 
increasing  values  of  time.  The 
opening  up  of  a  hole  in  the  back  of 
the  pulse  is  clear  from  both 
Figs.  18  and  19. 

Calculations  of  the  type  repre¬ 
sented  in  Figs.  17-19  become 
impractical  if  one  is  treating  a 


Fig.  17.  Saturation  of  on-axis  fluence 
due  to  strong  pulse  thermal 
blooming. 


-32- 


Fig.  18.  Three-dimensional  plot  of 
intensity  as  a  function  of 
time  and  radius  corresponding 
to  Figs.  16  and  17. 


Fig.  19.  Intensity  as  a  function  of 
radius  for  increasing  time 
in  pulse  corresponding  to 
Figs.  16  and  17. 


multipulse  beam.  The  determination 
of  nonisobaric  contributions  to  the 
density  is  greatly  simplified  by  the 
triangular  pulse  approximation,  in 
which  the  dependence  of  the  laser 
intensity  on  time  is  represented  as 
an  isosceles  triangle  with  base 
equal  to  2t^.  The  density  is 
required  only  at  time  t  =  T  ,  since 

tr 


the  laser  intensity  is  assumed  to 
vanish  for  t  =  0  and  t  -  2t  . 

r 

The  density  change  at  t  =  T  can 

r 

be  evaluated  analytically  in  terms 
of  a  finite  Fourier  series  represen¬ 
tation  of  the  laser  intensity.  The 
Fourier  transform  of  the  noniso- 
barically  induced  density  change  is 


(34) 


where  I  is  the  spatial  Fourier  transform  of  the  intensity.  The  corresponding 
density  changes  at  the  grid  points  are  given  by  the  discrete  Fourier  transform 
(DFT)  expression 

N 

p s*(jta,Uy)  =  (2»-2  £  pf  (=>  f )  W  (35) 

m,n=-N+ 1 
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where  the  basis  functions  are 
periodic  on  a  square  of  side  2 L. 

This  allows  for  a  buffer  region  that 
extends  an  additional  distance  L  in 
both  the  x  and  y  directions  from  the 
region  of  interest. 

Comparison  of  the  triangular 
pulse  approximation  and  detailed 
pulse  thermal-blooming  calculations 
for  Gaussian-shaped  pulses  in  time 
have  shown  good  agreement  between 
the  calculated  fluences  for  weak  or 
moderate  thermal  blooming.^- 


Fig.  20.  On-axis  fluence  as  a  function 
of  pulse  length,  as  calcu¬ 
lated  with  triangular  pulse 
approximation  (xfs)  and  by 
detailed  numerical  solution 
of  hydrodynamic  equations 
for  a  square  pulse  in  time 
(solid  curve).  The  tri¬ 
angular  pulse  approximation 
breaks  down  as  saturated- 
fluence  condition  sets  in 
at  T  =  1.5£  .  Erratic 
behavior  is  due  to  develop¬ 
ment  of  spikes  in  the 
intensity  pattern  as  a 
function  of  transverse 
position. 


Fig.  21.  Fluence  averaged  over 

minimum  area  containing  one 
half  of  total  beam  energy, 
as  a  function  of  pulse 
length.  Solid  curve  is 
detailed  calculation  for 
square  pulse,  xTs  represent 
triangular  pulse  approxi¬ 
mation. 

Figure  20  shows  the  on-axis  fluence 
calculated  for  the  previous  example 
with  the  triangular  pulse  approx¬ 
imation  (xfs)  and  the  detailed 
solution  of  Eqs.  (15)  for  square 
pulses  in  time  (solid  line) . 

Despite  the  difference  in  assumed 
pulse  shapes,  the  agreement  between 
the  two  types  of  calculation  is  very 
good  up  until  time  t  £2  50  ys,  which 
is  well  above  the  saturation  time 

t  =  38  ys  predicted  by  the  pertur- 

®  8 
bation  theory  of  Ulrich  and  Hayes 

9 

based  on  the  work  of  Aitken  et  al. 

Above  55  ys,  or  approximately  1.5t  , 

s 

the  beam  abruptly  develops  spikes 
in  its  transverse  spatial  dependence; 


-34 


this  clearly  signals  the  breakdown 
of  the  triangular  pulse  approxima¬ 
tion,  which  must  obviously  fail  when 
strong  saturation  behavior  sets  in. 

Figure  21  shows  the  fluence 
averaged  over  the  minimum  half¬ 
energy  area  (the  area  within  the 
one-half  peak  energy  contour)  cal¬ 
culated  with  the  triangular  pulse 
approximation  and  with  the  detailed 
solution  of  Eqs.  (15)  for  square 
pulses.  Both  calculations  increase 
initially,  reach  a  maximum,  and  then 
turn  over  with  increasing  time. 

This  is  in  part  due  to  the  increase 
of  the  area  within  the  one-half  peak 
energy  contour  with  time.  There  is, 
however,  no  point  in  believing  the 
triangular  pulse  approximation 
beyond  the  time  when  the  average 
fluence  curve  has  reached  a  maximum, 
which  also  coincides  with  the  onset 
of  erratic  behavior  in  the  on-axis 
fluence  (Fig.  20). 

The  perturbation  theory  alluded 
to  earlier  5  describes  the  on-axis 
fluence  saturation  for  a  beam  that 
is  initially  Gaussian  in  shape  and 
for  a  pulse  shape  that  is  square  in 
time.  In  this  theory,  the  expression 
for  the  on-axis  intensity  is 


'1 


nt)  =  i 


(s)i 


0 


t  <  t 

—  S 


(36) 

t>ts> 


where  ^q(s)  is  the  on-axis  intensity 
for  a  Gaussian  beam  propagating  in 
vacuum,  or 


I0U)  - 


_  -as 

IQ(0)  e 

D(z) 


(37) 


Here  a  is  the  absorption  coefficient 
and 


where  f  is  the  focal  distance  and  a 

is  the  radius  of  the  original 

Gaussian  beam.  The  saturation  time 

t  at  on- axial  position  z  is  given 
s 

by 


t  = 
s 


2N(y  -  1)  az2E  e  az 
_  LL - 


3'aa6D2(z)  T 


-1/3 


(39) 


where  N  is  the  refractivity ,  E 1  is 

the  pulse  energy,  and  is  the 

pulse  duration.  Since  the  fluence 

cannot  be  increased  for  pulses 

longer  than  t  ,  it  can  be  argued 
s 

that  nothing  is  accomplished  by 

making  the  pulse  longer  than  t  . 

s 

The  fluence  must  be  maximized 

instead  by  maximizing  the  product 

I Az)  t  or,  equivalently,  by 
u  s 

maximizing  (z) .  The  maximum 
allowable  value  of  1q(s)  at  point  z 
is  normally  determined  by  the  con¬ 
dition  that  it  not  exceed  the  break¬ 
down  intensity,  or 
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max  JQ(s)  =  JBD  • 


(43) 


(40) 


This  maximum  allowable  intensity  in 
turn  determines  a  critical  input 
pulse  energy  at  z  =  0  given  by 


E 


crit 


2,  t  N  az 
=  Tra  tsImD{z)e 


(41) 


where  Eq.  (37)  has  been  made  use  of, 

and  where  t  is  calculated  from 
s 

2N(y  -  1)  az2I 

t  =  - 7 - - 

s  3 a  D(z) 


-1/3 

.  (42) 


If  one  is  dealing  with  a  multipulse 
laser  with  pulse-repetition  fre¬ 
quency  V,  Eq.  (41)  can  be  used  to 
define  a  critical  input  power  with 


P  ..  =  VE  . 

crxt  crxt 

2  ,  N  az 

=  Tra  vtgIBDP(s)^ 

The  self-consistency  of  the 
triangular  pulse  approximation,  on 
the  other  hand,  prevents  the  on- axis 
intensity  from  ever  becoming 
negative,  but,  as  previously 
remarked,  the  triangular  pulse 
approximation  breaks  down  for  pulse 
energies  greater  than  the  value  that 
maximizes  the  space-averaged  target 
fluence.  For  this  pulse  energy,  the 
average  and  on-axis  fluences  should 
be  saturated,  and  further  increases 
in  pulse  energy  would  give  no  return. 
Figures  22  and  23  have  been  calcu¬ 
lated  with  the  data  on  which 


Fig.  22.  On-target  fluence  from  tri¬ 
angular  pulse  approximation 
averaged  over  area  contain¬ 
ing  (1  -  1/e)  fraction  of 
total  beam  energy.  Range 
=1.5  km,  JBD  =1.6 
x  106  w/cm^. 


Normalized  input  -  pulse  energy 

Fig.  23.  On-target  space-averaged 
fluence  and  intensity  as 
functions  of  input  pulse 
energy  for  triangular  pulse 
approximation.  Range 
=  2  km,  Pgp  =  3  x  10  W/cm^, 
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Figs.  16-21  are  based,  but  with  the 

following  differences:  the  ranges 

for  Figs.  22  and  23  are  1.5  km  and 

2.0  km  respectively;  the  values 

assigned  (somewhat  arbitrarily)  to 

6,2 

T  at  these  ranges  are  3  x  10  W/ cm 

bD  g  2 

and  1.5  x  10  W/cm  . 

Both  the  on-target  space-averaged 

fluence  and  intensity  (Fig.  23)  are 

plotted  as  functions  of  the  input 

pulse  energy  normalized  to  #crit 

given  in  Eq.  (42).  The  space 

averaging  is  over  the  area  contained 

within  the  1/e  energy  contour.  The 

indicated  maxima  of  the  average 

fluences  in  both  Figs.  22  and  23 

occur  at  an  input  pulse  energy  equal 

to  1.7 E  .  .  The  space-averaged 
crit 

fluence  curves  in  Figs.  22  and  23 
are  smoother  than  those  displayed 
in  Fig.  20  because  the  former  are 
averaged  over  larger  areas .  The 


scaling  implications  of  the  pertur¬ 
bation  theory  described  in 
Eqs .  (36)-(42)  are  apparently  valid 
for  the  triangular  pulse  approx¬ 
imation,  although  the  maximum  useful 
pulse  energy  predicted  by  the  latter 
is  about  50%  greater  than  that 
predicted  by  the  perturbation 
theory . 

In  summary:  the  triangular  pulse 
approximation  should  provide  reason¬ 
ably  accurate  fluence  results  for 
pulse  energies  up  to  the  values 
where  strong  thermal  blooming 
saturates  the  on-axis  fluence.  The 
breakdown  of  the  approximation  will 
be  indicated  by  the  development  of 
spikes  in  the  transverse  spatial 
dependence  of  the  beam  intensity  as 
well  as  by  a  sharp  falloff  in  the 
fluence  averaged  over  some  area  as 
a  function  of  pulse  energy. 


11.  Multipulse  Thermal  Blooming  in  the 
Triangular  Pulse  Approximation 


The  propagation  of  a  given  pulse 
in  a  train  is  influenced  by  both  the 
nonisobaric  density  changes 
discussed  in  the  previous  section 
and  by  the  isobaric  density  changes 
due  to  heating  by  previous  pulses 
in  the  train.  But  can  the  self¬ 
blooming  and  multipulse  blooming 
effects  be  treated  independently? 


If  so,  the  results  and  discussion 
of  the  previous  section  suggest 
that,  as  time-averaged  laser  power 
is  increased  by  lengthening  the 
duration  of  the  constituent  pulses 
in  the  train,  the  time-averaged 
intensity  on  target  should  saturate 
at  a  value  that  is  predictable  from 
the  saturation  fluence  for  a  single 
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pulse.  If  <I>  represents  the  time- 
averaged  intensity,  the  maximum 
achievable  value  of  for  a  given 
pulse-repetition  rate  should  be 
expressible  as 


<J>  =  VF  _  , 

max  sat  5 


(44) 


where  F gat  is  the  single-pulse 
saturation  fluence. 

In  order  to  test  the  hypothesis 
of  the  independence  of  self  and 
multipulse  blooming,  a  set  of  cal¬ 
culations  has  been  carried  out  with 
the  following  set  of  parameters: 
Start  beam  shape  Gaussian,  truncated 


at  lie  radius 
Range,  R 
Focal  length/ 
range,  F/R 
Wavelength,  X 
Absorption 

coefficient,  a 
Aperture  diameter, 
2 a  (Gaussian  at 
l/e2) 

Wind  velocity,  v q 

Pulse-repetition 

rate,  V 

Maximum  pulse 

intensity  at 

receiver,  I 

max 

Overlap  number, 

NQ  =  2cco/ Vq 

Figure  24  shows 
averaged  single-puj 
for  V  =  33-1/3  s-1 


2.5  km 

1.0  and  1.2 
10.6  ]im 

0.25  km  ^ 

21.2  cm 
10  m/s 

33-1/3  and  50  s 

4.9  MW/ cm2 

1.0,  1.5 

the  space- 
se  intensity  I 
and  Nq  =  1,  with 


F/R  -  1.0  and  1.2,  calculated  as  a 
function  of  input  time-averaged  power 
<P>  =  \)Ep.  The  curves  have  been  cal¬ 
culated  with  and  without  the  effects 
of  pulse  self-blooming.  The  curve 
without  self-blooming  for  F/R  =1.2 
rises  slightly  with  input  power 
because  of  a  very  slight  amount  of 
pulse  overlap.  It  is  clear  from 
Fig.  24,  in  any  case,  that  thermal 
blooming  is  due  almost  entirely  to 
self-blooming  effects.  The  corres¬ 
ponding  curves  for  space-  and  time- 
averaged  target  intensities  <J>  are 
displayed  in  Fig.  25,  where 

<J>  =  Jtv  .  (44) 

It  is  seen  that  <I>  with  self-blooming 
rises  initially,  reaches  a  peak,  and 
then  falls.  From  the  analysis  of  the 
previous  section,  we  interpret  the 
peak  values  of  <J>  as  the  saturated 
values . 

Figure  26  shows  I  as  a  function 
of  <P>  for  V  =  50  s-1  and  NQ  =  1.5, 
with  F/R  =  1.0  and  1.2.  Above  <P>  = 
=0.5  MW,  and  an  enhancement  effect 
sets  in  that  is  greater  in  the  case 
of  the  defocused  beam.  The  corres¬ 
ponding  curves  for  space-  and  time- 
averaged  target  intensities  are  shown 
in  Fig.  27. 

A  comparison  of  Figs.  25  and  27 
is  summarized  in  Table  8.  It  is  seen 
that  at  V  =  50  s_1  the  power  <P> 

sat 

at  which  saturation  of  <J>  occurs  is 


-38- 


target  intensity  —  kW/ 


Fig.  25. 


Space-  and  time-averaged  intensity  on  target  as  a  function  of  time- 

averaged  power  at  transmitter:  v  =  33-1/3  s-1  Nn  =  1.  (a)  F/R  =  1  0 

w/n  =  io  ’  0  v  -L*u* 
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Space-averaged  target  intensity  —  kW/ 


8 


Fig.  26.  Space-averaged  intensity  on  target  as  a  function  of  time-averaged 
power  at  transmitter:  V  =  50  s~l ,  No  =  1.5.  (a)  F/R  -  1.0. 

(b)  F/R  =  1.2. 
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Space-  and  time-averaged  target  intensity  —  kW/< 


Time-averaged  transmitter  power  —  MW 

Fig.  27.  Space-  and  time-averaged  intensity  on  target  as  a  function  of  time- 
averaged  power  at  transmitter:  v  =  50  s~l  Nn  =  1  5 
(a)  F/R  =  1.0.  (b)  F/R  =  1.2,  ’  ° 
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Table  8.  Saturation  of  time-  and 
space-averaged  target 
intensity  due  to  self- 
blooming. 


V 

(s"1) 

F/R 

<P>  - 

sat 

(MW) 

<I>  - 
sat 

(kW/cm^) 

33  1/3 

1.0 

1.0 

2.4 

33  1/3 

1.2 

1.2 

2.5 

50 

1.0 

1.75 

2.7 

50 

1.2 

2.0 

3.5 

higher  for 

both 

values  of  F/R  than 

it  is  at  V 

=  33- 

-1/3  s  1. 

The  cor- 

responding  saturation  intensity 

values  <I>  are  also  greater  at 

-1Sat  -1 
V  =  50  s  than  at  V  =  33-1/3  s 

If  effects  of  self-blooming  are  not 

included,  on  the  other  hand,  values 

of  <J>  are  always  greater  at  a  given 


value  of  <P>  in  the  case  of 
V  =  33-1/3  s"1. 

Unfortunately,  we  have  no  guide 
to  the  accuracy  of  the  triangular 
pulse  approximation  in  the  overlap 
case  as  we  do  in  the  nonoverlap  case. 
But  the  above  results  strongly  sug¬ 
gest  that  the  contributions  of  iso- 
baric  and  nonisobaric  density  changes 
to  thermal  blooming  of  multipulse 
beams  are  interrelated,  and  that 
time-averaged  saturation  intensities 
based  on  single  saturation  fluences 
may  not  be  applicable  for  overlap 
numbers  somewhat  above  1.  In  fact, 
overlapping  isobaric  density  patterns 
may  in  certain  situations  actually 
override  the  effects  of  single-pulse 
nonisobaric  density  changes. 
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Appendix  A:  Adaptive  Lens  Transformation 


One  key  to  the  successful  implementation  of  a  laser-propagation  code 
is  finding  a  coordinate  transformation  that  keeps  the  laser  beam  away  from 
the  calculational  mesh  boundary  and  at  the  same  time  prevents  the  beam  from 
contracting  to  an  unreasonably  small  fraction  of  the  total  mesh  area  at  the 
focus.  If  one  is  solving  the  Fresnel  equation  by  the  finite  Fourier  trans¬ 
form  method,  one  may  alternatively  view  the  problem  in  terms  of  comple¬ 
mentarity:  one  wishes  to  find  a  transformation  that  simultaneously  keeps 

the  beam  intensity  small  on  the  mesh  boundaries  in  configuration  space  and 
keeps  the  Fourier  spectrum  small  on  the  mesh  boundaries  in  k-space.  If 
these  two  conditions  are  met,  one  knows  from  sampling  theory  that  the 
numerical  solution  is  highly  accurate. 

The  Four-D  code  uses  an  automated  procedure  that  is  designed  to  keep 
the  intensity  centroid  at  the  center  of  the  mesh  and  the  intensity-weighted 
r.m.s,  values  of  x  and  y  constant  with  propagation  distance  z.  These  con¬ 
ditions  can  be  written 


> 

I 


0, 


(Ala) 


9 

9a 


<xi  ~  <xi>> 


2 

> 

I 


=  0  , 


i  =  1,  2  ,  (Alb) 


x±  =  x,  x2  =  y  , 


where 


<u> i  = 


./■ 


dx  dy  lCx3ylu 


/ 


dx  dy  T(xsy) 


(Ale) 


Hereafter,  all  averages  will  be  assumed  to  be  intensity-weighted,  and  the 
subscript  I  will  be  dropped. 

Conditions  (Al)  also  apply  to  the  adaptive  coordinate  transformation 
A1 

of  Bradley  and  Hermann,  which  differs  from  the  one  employed  in  the  Four-D 


Al.  L.  C.  Bradley  and  J.  Hermann,  "Change  of  Reference  Wavefront," 

Massachusetts  Institute  of  Technology,  Lincoln  Laboratory,  Lexington, 
Mass.,  unpublished  internal  report. 
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code  only  in  that  it  is  preceded  by  a  transformation  to  the  coordinates  of 
an  arbitrary  Gaussian  beam  propagating  in  vacuum.  It  should  be  evident,  in 
any  case,  that  such  adaptive  transformations  are  restricted  to  steady-state 
problems,  since  for  time-dependent  problems  no  single  transformation  will 
apply  to  all  time  values.  To  solve  time-dependent  problems  one  must  employ 
a  Talanov  transformation  that  is  optimized  to  all  time  values.  This 
optimization  is  accomplished  by  a  combination  of  trial  and  error  and 
intuition. 

The  splitting  algorithm  employed  in  the  Four-D  code  can  be  written 
formally  as 


t?n+l  /  jJsz  n2\  /  -\  /  -£As  2\ 

9  -  exp  (-  jjp  VjJ  exp  (-  X)  exp  (-  ^  vj  t„.  <“) 

X  =  k2(n2  -  1)  , 

where  the  middle  exponential  on  the  right-hand  side  of  Eq.  (A2)  contains  the 
changes  in  phase  resulting  from  hydrodynamic  changes  in  density,  turbulence, 
etc.  Immediately  after  this  step  in  the  calculation,  a  quadratic  reference 
phase  front  is  determined  and  is  removed  from  '&  by  means  of  a  Talanov 
transformation  and  a  deflection  of  the  beam  coordinates.  These  operations 
are  carried  out  as  part  of  the  vacuum  propagation  step.  During  vacuum 
propagation  the  solution  is  advanced  by  solving 

Hk  H  =  Vi  8  •  (A3) 


(A4a) 

i  =  1,  2  (A4b) 

where  ~P  is  the  beam  power  given  by 

&x1  dx2  |^(a?1,a?2)  I  .  (A5) 


Equation  (Al)  can  be  written 


<>-?/ 


<Xi>  =  p  I  ^2  x±\$(xi>x2>b)\  9 


<x 


2  1  /  2  2 
^  “  p  I  I  9 
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By  differentiating  Eqs.  (A4a)  and  (A4b)  with  respect  to  s  and  making 
use  of  the  Fresnel  equation  (A3),  one  obtains  the  following  relations: 


9  ^  2^_  2_ 

3s  <xi>  kP 


^  (3^2  , ^2  > ^0  I  x  3^  d  (^2.  *  x2. 9  9 


-5—  <XJ 
dS  % 


t>--fe/*lfc2  „(W> 


2  3^/ 

3x.  4>(^i»«2’  ^ 


where  the  phase  (J>(rc^,  a^*  3)  i-s  defined  by 

4> »£C2 * s)  =  Im[ln  S^x-jL>x29Z^^  ' 


In  k-space  one  can  similarly  derive 


(A6a) 

(A6b) 


(A6c) 


-5—  <x  ■>  = 
3s  1 


2  <Ki> 


3  .  2  _  2  T 

3s  Xi  >  kP  ^ 


kP  JJ  dlcl  dK2  Kil^KrK2’2)|  k 

^ ^  d^2  di<2  kn  ,  i<2  j^)  I 


(A7a) 


X  3k“  t(K1»<2’s^  »  (A7b) 

i 

'■w 

where  ^(k^j^jS)  is  the  Fourier  transform  of  & (a?^,aj^,s)  ,  and 

iKk15K2,s)  =  Im[ln  (ic^i^s)]  .  (A8) 

We  now  wish  to  determine  a  phase  front  that  will  preserve  the  following 
conditions: 


Oc  •>  =  0  ,  (A9a) 

|i-<(iC  -<xi»2>=  0  ,  i  =  1,  2.  (A9b) 
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Equation  (A9b)  is  equivalent  to 


3  2  „ 

■5—  <x  .  >  -  2<x  •> 

03  1  1 


■5—  <x.>  =  0. 

03  1 


From  Eqs.  (A9)  and  (A6)  one  obtains 


Thus  the  reference  phase  front  must  satisfy 


Let  us  define  a  new  phase  variable: 


(A9c) 

(AlOa) 

(AlOb) 

(Alla) 

(Allb) 


^(xrxvzn^  =  VVVW 


+  J  [a^  (aN  -  <x^>)2  +  -  <xx>)]  , 


(A12) 


where  <J>q( x\>x2’Zn\%?  represents  the  phase  at  Sn-^  before  the 

vacuum  propagation  operator  has  been  applied,  and  where  0^  and  are 
determined  so  as  to  make  conditions  (Alla)  and  (Allb)  hold  for  the  phase 
front  (x1>x2’zn+k*  *  From  Ec5,  (Alla)  we  obtain 


\ST  -w)  "  2a<  ‘ 


<#.»+$. 

t-  1. 


+  V^7  ♦0(arl»*2»  Wy 
%  • 


(A13) 
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or 


From  (Allb)  we  obtain 


<^>) 


dx  . 

1 


<f>*  («i 


,x2  ,sn+V 


or 


+  3  .<x .  -  <x .»  =  0 

i  i  i 


i  =  1,  2  .  (A14) 


Equations  (A13)  and  (A14) ,  which  determine  the  desired  reference 

phase-front  parameters  (A12),  can  be  shown  to  be  completely  equivalent  to 

A1 

the  relations  used  by  Bradley  and  Hermann. 

If  the  optimal  phase  front  (J>'  is  now  substituted  for  the  original 
phase  front  <J>Q  at  sn+ig,  the  phase  increment 


2 

<(>q  -  <t>'  =  -  ^  ~  ~  (^15) 

i- 1 


be  compensated  for  xn  some  way  in  order  to  preserve  the  original  field. 
The  quadratic  contribution  in  (A15)  is  compensated  for  by  a  generalized 
Talanov  transformation,  which  involves  a  rescaling  of  <?,the  mesh,  and  As, 
according  to 


$(,x ,  y ,  s) 


1 


x  S 


x 

As  ’ 


exp 


(A16a) 
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where 


£(Kx’Ky’s^  =  ^Kx*Ky*8  “  As)  exp 


(A16b) 


z 


1 


_ z 

1  -  (A s/s  ) 

X 


(A16c) 


_ z 

1  -  (EzJlT) 


(A16d) 


The  generalized  focal  lengths  and  z are  determined  by  combining  the 
reciprocals  of  the  current  focal  lengths. 


=  2 Oj/fc  , 

(A17a) 

=  2a2/k  , 

(A17b) 

wi-th  those  remaining  from  previous  propagation  steps  (see  argument  of  expo¬ 
nential  in  Eq,  (A16a)). 

The  linear  term  in  (A15)  corresponds  to  solving  Eq.  (A3)  in  a  coor¬ 
dinate  system  that  has  been  rotated  in  x-y-z  space.  If  this  rotation  is 
assumed  to  be  small,  it  can  be  represented  by  a  net  deflection  in  the  x  and 
y  coordinates  given  by 


6a;  =  -(31/fc)g  , 


(A18a) 


6y  =  -(3 2fk)z  . 


(A18b) 


The  contribution  must  also  be  added  to  <pQ  before  the  vacuum 

propagation  calculation,  but  this  operation  may  correspond  to  a  translation 
of  the  Fourier  transform  )  by  a  nonintegral  number  of  steps  on  the 
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k- space  mesh.  In  order  to  avoid  this,  3-,/Ak  and  39/Ak  are  both  rounded 

j_  cc  £  y 

off  to  the  nearest  integer,  and  k  )  is  then  translated  on  its  mesh  in 

x  y 

the  x  and  y  directions  by  the  corresponding  number  of  steps. 

The  numerical  implementation  of  Eqs,  (A13)  and  (A14)  requires  the 
following  computations,  where  j  and  k  represent  the  numerical  coordinates 
of  the  mesh  points: 


<a;>  E  J  Xj\Sjk}  * 
3  ,k 


(ix  <x>)  y  E  ^  i 

3  ,k 


x 


(A19) 


tk  j  ^ k 

3  = 


x  \  i  ~  1 2 

0\ 


’  =  (2Aa;ffrl  Im  J  %-k  ~  %-l,k)(4k  +  4-1^  , 

0,k 


{(x  -  <x>)  -  ( 2bxE) 


-1 


X  Im  S  x3  (<?J k  “  *j-l ,fc)(<gj7c  +  *3-1  ,k)  <x>  (H) 

3,k 
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The  computations  involving  the  variable  y  are  carried  out  in  an 
analogous  manner.  In  the  calculation  of  the  average  phase  derivative,  the 
phase  derivative  is  monitored  at  each  point  and  limited  in  magnitude  to  a 
fraction  of  ir.  This  prevents  rapid  phase  fluctuations  near  the  mesh 
boundary,  where  intensities  may  be  weak,  from  contributing  disproportionately 
to  the  average. 
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Appendix  B:  An  Adaptive  Algorithm  for  Selecting  the 

Axial  Space  Increment 

It  is  desirable  to  have  the  code  select  the  next  axial  space  increment 
A z  at  a  given  axial  position  on  the  basis  of  requirements  for  numerical 
accuracy  in  the  solution  of  the  wave  equation.  The  numerical  acquracy  of 
the  vacuum  propagators  in  the  symmetrically  split  solution  operator, 

*>n+ 1  /  ihz  n.  2\  /  ihz  -\  (  ihz  „  2\  ^n 

S  =  exp  I  -  Vj_  J  exp  I-  2^-  xl  exp  I-  V_|l<f,(Bl) 

is  independent  of  A z  if  the  solution  is  based  on  a  discrete  Fourier  trans¬ 
form.  The  imposition  of  the  phase  front, 

Af-#  ,  (B2) 

at  z  =  which  is  equivalent  to  passing  the  beam  through  a  lens,  will 

make  the  solution  meaningless  if  any  of  the  transverse  zone-to-zone  phase 
differences  violate 

|6^(A4>)|  <  frr  ,  (B3) 

|  4  /it  , 

0  <  f  <  1  . 

It  will  always  be  necessary  then  to  restrict  the  value  of  As  so  that  con¬ 
ditions  (B3)  are  met.  While  violating  conditions  (B3)  destroys  the  numerical 

integrity  of  the  solution,  satisfying  them  does  not  completely  guarantee 

2 

accuracy,  since  errors  can  also  result  from  the  noncommutation  of  Vj_  and 
X,  and  from  upgrading  x  to°  infrequently.  These  errors  must  be  controlled 
externally  by  inputting  a  maximum  allowable  value  of  As. 

In  practice,  part  of  the  effect  of  the  phase  front  (B2)  is  removed  by 
the  adaptive  lens  transformation.  It  would  therefore  be  too  restrictive  to 
limit  Ajs  on  the  basis  of  conditions  (B3),  As  an  alternative  one  can  restrict 
the  value  of  A z  so  as  to  control  transverse  gradients  in  the  phase  variable 
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S  ai^Xi  ~  <Xi>-?2  ”  2  3p  As npl  ’ 


which  is  that  part  of  the  nonlinear  phase  front  at  a  tl  that  cannot  be 

n+% 

removed  by  the  adaptive  lens  transformation.  The  next  spacial  increment 

77.~4*1  ,  yi 

As  is  then  chosen  in  terms  of  the  current  value  As  by  means  of  the 
relation 


, _ I**  71 _ , _  /tj  r  \ 

max  110  I  ~  -  0  .7  I,  10  ~  7  -0.-1}  # 

v  1  jk1’  1  j+l,fe  >f»j 

A:—  max 


The  arguments  of  the  maximum  function  in  the  denominator  of  expression  (B5) 
are  restricted  to  those  mesh  points  where  the  intensity  is  greater  than  a 
certain  fraction  of  the  maximum  intensity.  The  final  value  of  Asn+\ 
however,  must  satisfy  the  additional  constraints: 

0.8Asn  <  Aan+1  <1.2  A zn  ,  (B6) 


A z  .  <  A zn+^  <  As 

mxn  —  max 


9 


+ 


k\_  (^x  -  <x>)2y .+  ({y  -  <i/>)2)] 


(B7) 

(B8) 


with  (B6)  taking  precedence  over  (B5),  (B7)  over  (B6) ,  and  (B8)  over  (B7) . 
In  Eq.  (B8),  fn  *»  0.005  is  an  input  fraction  and  s„  =  min(|s  |,  \z  |). 
Condition  (B8)  is  designed  to  reduce  As  near  a  focus,  where  the  geometric- 
optics  scaling  of  the  mesh  by  the  Talanov  transformation  may  result  in  an 
excessive  shrinkage  of  the  mesh.  By  updating  the  Talanov  transformation 
sufficiently  often,  one  can  usually  avoid  a  geometric-optics  catastrophe. 

The  adaptive  s-step  algorithm  just  described  adds  greatly  to  the  con¬ 
venience  of  running  problems;  it  often  improves  problem  running  time,  and 
avoids  large  nonlinear  phase  changes  that  can  invalidate  the  calculation. 

It  should  not,  however,  be  regarded  as  a  panacea.  For  sufficiently  high 
beam  power  and  strong  enough  thermal  blooming,  the  criteria  (B5)-(B8)  can 
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be  satisfied  and  yet  the  problem  still  goes  bad.  In  such  cases,  large  non¬ 
quadratic  transverse  zone-to-zone  phase  differences  can  accumulate  over  many 
2-steps. 
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Appendix  C:  Treatment  of  Multiline  Absorption 

The  treatment  of  multiline  absorption  in  the  Four-D  code  follows  the 
Cl 

method  of  Hogge.  The  basic  assumptions  are  that  all  lines  operate  with 
the  same  transverse  mode  structure  in  the  laser  and  that  the  line  frequen¬ 
cies  are  near  enough  to  each  other  so  that  the  field  for  each  line  will  be 
affected  in  the  same  way  by  the  atmospheric  density  distribution.  Thus  at 
each  position  3,  for  the  ith  line  <?.(a)  cc  S{z)  ,  and  the  field  will  be  com- 

Is 

pletely  characterized  by  the  fractions  f .  (s)  of  the  total  power  P(z )  that 

is 

are  found  in  each  line. 

At  z  =  0  one  has 


P.(0)  =  /\( 0)  P(0)  , 


and  at  position  z 


P^z)  =  Pt{ 0)  a  aiz  , 


F(z)  =  2  FAz)  =  p(°>  6 


-a .  z 


fAz)  = 


fA °)  e 


■a  .z 


2  ^(0)  e 


-a  ,z 

i 


Thus, 


jAs)  ~  fA3'*  l^)> 


and  the  energy-deposition  rate  per  unit  volume  is  given  by 


(Cl) 


(C2) 


CC31 


Cl.  C.  B.  Hogge,  "A  Comparison  of  Several  High  Energy  Laser  Systems  with 
Emphasis  on  Propagation  Aspects,"  in  Laser  Digest,  AFWL-TR-75-140 
(May  1975). 
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(C4) 


2aiJi(s)  =  J(s)  5 

i  i 

in  which 

=  (C5) 

i 

can  be  interpreted  as  an  average  cross  section  throughout  the  calculation. 

Tables  Cl  and  C2  show  values  for  a(z)  as  a  function  of  z  for  the  DF 
line  data  found  in  Ref.  C2.  Clearly,  for  DF  the  effect  of  including  all 
line  absorption  details  leads  to  a  very  small  correction  even  at  10  km. 


Table  Cl.  Line-by-line  absorption-coefficient  data. 


Line  ID 

Line  frequency 
(cm~l) 

Fraction  of 
2  =  0 

total  power 

Z  —  10  km 

Absorption 
coefficient 
(km-1  ) 

4-3,7 

— 

0.01040 

0.01006 

0.06000 

3-2,10 

2496.77 

.00590 

.00636 

.04920 

4-3,6 

— 

.02130 

.02060 

.06000 

3-2,9 

2521.81 

.01330 

.01632 

.03620 

4-3,5 

— 

.01040 

.01006 

.06000 

3-2,8 

2546.42 

.04750 

.05401 

.04380 

3-2,7 

2570.51 

.06380 

.06108 

.06100 

2-1,10 

2580.10 

.00910 

.00813 

.06790 

3-2,6 

2594.25 

.08970 

.11944 

.02800 

2-1,9 

2605.80 

.03180 

.03786 

.03920 

3-2,5 

2617.44 

.05630 

.07960 

.02200 

2-1,8 

2631.06 

.08450 

.10460 

.03530 

2-1,7 

2655.85 

.09040 

.05918 

.09900 

2-1,6 

2680.17 

.13040 

.12864 

.05800 

1-0,9 

2691.61 

.03230 

.02557 

.08000 

2-1,5 

2703.99 

.04000 

.05300 

.02850 

1-0,8 

2717.54 

.06440 

.03629 

.11400 

1-0,7 

2743.00 

.08740 

.08434 

.06020 

1-0,6 

2767.97 

.08370 

.06178 

.08700 

1-0,5 

2792.43 

0.02740 

0.02310 

0.07370 

C2.  R.  K.  Long,  F.  S.  Mills,  and  G.  L.  Trusty,  Calculated  Absorption- 

Coefficients  for  DF  Laser  Frequencies,  Ohio  State  University  Electro- 
Science  Laboratory,  Rept.  RADC-TR-73-389 . 
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Table  C2.  Mean  absorption  coefficient  as  a  function  of  distance. 


Mean 

Remaining 

Propagation 

absorption 

power 

distance 

coefficient 

fraction 

(km) 

(km“l) 

in  beam 

0 

0.06001 

1.00000 

1 

0.05931 

0.94209 

2 

0.05862 

0.88815 

3 

0.05793 

0.83787 

4 

0.05726 

0.79097 

5 

0.05660 

0.74720 

6 

0.05594 

0.70632 

7 

0.05530 

0.66811 

8 

0.05467 

0.63236 

9 

0.05404 

0.59891 

10 

0.05343 

0.56758 

I 
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Appendix  D:  Characterization  of 
Nondiffraction-Limited  Beams 

In  the  absence  of  detailed  a  priori  information  regarding  the  exact  mode 
content  of  a  beam,  several  models  of  nondiffraction-limited  beam  behavior  can 
be  applied  with  the  Four-D  code. 

The  simplest  of  these,  which  requires  no  special  coding,  is  wavelength 
scaling,  wherein  the  laser  wavelength  is  multiplied  by  a  number  equal  to  the 
beam  quality  factor.  Wavelength  scaling  gives  the  correct  vacuum  peak 
intensity,  although  it  may  incorrectly  represent  the  vacuum  focal-spot  size. 

It  represents,  in  any  case,  a  prescription  whose  accuracy  needs  to  be  evaluated 
ad  hoc  for  each  specific  application.  While  it  has  been  useful  in  a  variety 
of  applications,  it  does  not  properly  account  for  discrepancies  between  cal¬ 
culations  and  stagnation sb looming  experiments  in  vertical  absorption  cells. 

Agreement  between  measured  data  and  calculations  for  these  experiments 
is  improved,  on  the  other  hand,  by  adding  spherical  aberration  to  the  initial 
beam  in  such  a  way  that  the  vacuum  focal-spot  size  is  correctly  reproduced 
(see  Fig.  Dl) .  The  spherical  aberration  contribution  to  the  initial  phase  can 
be  represented  as 


(Dl) 


where  A  represents  the  number  of  waves  of  aberration  at  radius  0  . 

D2 

A  third  model  of  nondiffraction-limited  behavior  is  due  to  Hogge  et  dl . 
This  model  is  based  on  the  assumption  that  the  initial  beam  can  be  represented 
by 

G(x,y,0)  =  Sq(x,u)  ,  (D2) 

where  the  phase  aberration  $(x9y)  is  a  Gaussian  random  variable,  arising  from 
laser-medium  inhomogeneities,  mirror  imperfections,  etc.  If  it  is  assumed 
that  the  correlation  function  for  phase  fluctuations  is 


Dl.  J.  A.  Fleck,  Jr. ,  J .  R.  Morris,  and  M.  D.  Feit,  Time- Dependent 
Propagation  of  High  Energy  Laser  Beams  Through  the  Atmosphere , 

Lawrence  Livermore  Laboratory,  Kept.  UCRL-51826  (1975). 

D2.  C.  B.  Hogge,  R.  R.  Butts,  and  M.  Burlakoff,  Appl .  Opt.  13,  1065  (1974). 


Relative  intensity 


C.(r)  =  a 

$ 


2 

<f> 


exp 


(D3) 


where  =  x2  +  y2 ,  o2^  is  the  variance  of  <j),  and  l ^ 
length,  then  the  spectrum  of  the  phase  fluctuations 


is  the  phase  coherence 
is  given  by 


W  - 


2  72 

2t r 


exp  I  - 


i2  k2) 


(D4) 


Fig.  Dl.  Intensity  on  target  after  passing  through  stagnation  zone.  Com¬ 
parison  between  experiment  and  calculation  with  and  without  spherical 
aberration.  (Data  from  P.  J.  Berger,  F.  G.  Gebhardt,  and  D.  Smith, 
Thermal  Blooming  Due  to  a  Stagnation  Zone  in  a  Slewed  Beam ,  United 
Aircraft  Research  Laboratory,  East  Hartford,  Conn.,  Rept.  N921724-12 
(1974).) 
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Using  the  method  of  phase  screens,  one  can  obtain  the  Fourier  transform  of 
the  phase  §(,x}y)  to  be  used  in  Eq.  (Dl)  in  the  following  form: 


ta'  (kp  +  id'  (k^ 

V  VT  ) 


(D5) 


where  a'  and  a"  are  Gaussian  random  variables  with  variance  1,  and  where 


a'<V  =  a'(’V  ’ 

(D6) 

a"<V  ’  • 

Equations  (D5)  and  (D6)  were  originally  included  in  the  Four-D  code  for 
simulating  turbulence,  The  phase-screen  model  of  nondiffraction-limited 
beams  utilizes  the  same  subroutines. 

The  following  parameters  were  the  basis  of  an  example  for  comparing  the 


difference  between  the  wavelength-scaling  and  the  phase-screen  models  of  non- 

diffraction-limited  behavior: 

Beam  shape 

Gaussian 

Aperture  size,  2 a 

80  cm 

Range 

2. 5  km 

Absorption  coefficient,  a 

0.07  km 

Transverse  wind  speed,  V 

10  m/s 

Focal  distance,  / 

Beam  is  2 *  diffraction  limited 

2.5  km 

Wavelength,  X 

5.7  ym 

Scaled  wavelength,  X 

11.4  ym 

Phase  correlation  length,  l ^ 

5 . 0  cm 

Phase  standard  deviation,  0 ^ 

Number  of  phase-screen  calculations 

1.177  (rad) 

for  ensemble  average 

10 
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Table  D1  gives  a  comparison  of  peak  intensities  in  the  focal  plane  for  propa¬ 
gation  in  vacuum  and  air.  Also  included  are  results  for  a  uniformly  illumina¬ 
ted  aperture  of  radius  cZq  =  2 a,  which  can  likewise  be  used  as  a  model  of  a  2x 
diffraction-limited  beam. 


Table  Dl.  Comparison  of  peak  intensities  in  the  focal  plane.  Beam  propagates 

in  vacuum  (linear)  and  air  (nonlinear).  Case:  Gaussian  diffraction- 
limited  beam,  a  uniformly  illuminated  aperture  (top  hat),  Gaussian 
wavelength  scaled  (2*  diffraction  limited) ,  and  a  Gaussian  beam 
with  a  phase  screen  adjusted  to  2x  diffraction  limited. 


Model 

Peak  intensity  in 
focal  plane  (kW/cm^) 

Linear 

Gaussian,  A 

63.8 

Gaussian,  2A 

15.9 

Top  hat 

17.4 

Phase  screen 

17.7 

Nonlinear 

Gaussian,  2X 

5.94 

Top  hat 

5.47 

Phase  screen 

3.03 

The  nondiffraction-limited  beams  all  give  roughly  one  quarter  of  the 
focal-plane  intensity  of  the  diffraction-limited  beam  when  propagated  in 
vacuum.  For  propagation  in  a  real  absorbing  atmosphere  both  the  scaled  wave¬ 
length  and  the  top-hat  beam  calculations  result  in  twice  the  peak  intensity 
of  the  phase-screen  model  calculation,  which  is  based  on  an  ensemble  average 
taken  over  10  independent  phase  screens.  Figures  D2  and  D3  show  respectively 
the  beam  intensity  as  a  function  of  position  along  the  ar-axis  and  along  a  line 
parallel  to  the  y- axis  passing  through  the  point  of  maximum  intensity  along  the 
a>axis  for  the  wavelength-scaled  beam.  Figures  D4  and  D5  show  the  same  plots 
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Fig.  D2.  Wavelength  scaling  for 

2x  diffraction-limited  beam. 
Intensity  on  target  as  a 
function  of  x  along  rc-axis. 


Fig.  D4.  Phase-screen  model  of  2x 
diffraction-limited  beam 
averaged  over  10  independ¬ 
ent  realizations.  Inten¬ 
sity  on  target  as  a 
function  of  x  along  rs-axis. 


Fig.  D3.  Wavelength  scaling  for  2x 
diffraction-limited  beam. 
Intensity  on  target  as  a 
function  of  y  along  a  line 
parallel  to  y-axis  and 
passing  through  point  of 
maximum  intensity  along  the 
re-axis. 

for  the  ensemble  averaged  phase- 
screen  calculation.  In  the  case  of 
the  wavelength-scaled  calculation, 
thermal  blooming  leads  primarily  to 
a  broadening  of  the  beam.  In  the 
phase-screen  calculation,  thermal 
blooming  is  accompanied  by  con¬ 
siderable  scattering  of  energy  to  the 
far  reaches  of  the  mesh. 

From  these  calculations  it  can 
be  concluded  that  different  models 
of  nondiffraction-limited  beam 
behavior  can  lead  to  qualitatively 
as  well  as  quantitatively  different 
thermal-blooming  behavior.  Which 
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Intensity  —  kW/cm 


model  is  best  must  be  determined  for 
each  specific  situation.  In  the  last 
analysis  there  is  no  substitute  for  an 
accurate  experimental  characterization 
of  the  beam  for  each  specific  laser. 

Fig.  D5.  Phase-screen  model  of  2x 
diffraction-limited  beam 
averaged  over  10  independ¬ 
ent  realizations.  Inten¬ 
sity  on  target  as  a 
function  of  y  along  a 
line  parallel  to  z/~axis 
and  passing  through  point 
of  maximum  intensity 
along  the  x-axis. 
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